{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Issue 12](https://github.com/alkaline-ml/pyramid/issues/12)\n",
    "\n",
    "The purpose of this notebook is to debug issue 12\n",
    "\n",
    "Here is the Pyramid version info:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pyramid version: '1.0.0-dev'\n"
     ]
    }
   ],
   "source": [
    "import pmdarima as pm\n",
    "print(\"Pyramid version: %r\" % pm.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the data in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>occupancy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2017-03-01 00:02:01</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2017-03-01 00:04:01</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2017-03-01 00:06:01</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2017-03-01 00:08:01</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2017-03-01 00:10:01</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  time  occupancy\n",
       "0  2017-03-01 00:02:01          2\n",
       "1  2017-03-01 00:04:01          3\n",
       "2  2017-03-01 00:06:01          2\n",
       "3  2017-03-01 00:08:01          1\n",
       "4  2017-03-01 00:10:01          4"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('dummy_data.csv')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determining the periodicity\n",
    "\n",
    "Our issue filer claims that the periodicity of the data is every 2 minutes:\n",
    "\n",
    "> I have a seasonal time-series data set which comprises of 2 months of data. Now if I mention the frequency as \"2 min\" or \"365x24x30\" then it tries to divide the dataset with that frequency. But my data set is just for 2 months and obviously, it has way lesser values than the freq.\n",
    "\n",
    "Therefore, rather than use `365 * 24 * 30 = 262800`, we can probably use the number of days in the two months (say, 60) instead of the 365: `60 * 24 * 30 = 43200`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43200"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_days = 60\n",
    "n_hours = 24\n",
    "n_minutes = 30  # every other minute\n",
    "\n",
    "# determine m:\n",
    "m = n_days * n_hours * n_minutes\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, `m` is greater than the number of samples in the data, which is suspect... *unless there is LOTS of data, it's unlikely we have a bi-minutely seasonal period...*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_periods: 43200\n",
      "n_samples: 1610\n"
     ]
    }
   ],
   "source": [
    "print(\"n_periods: %i\" % m)\n",
    "print(\"n_samples: %i\" % data.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On inspection of the data, it's clear that there is a 2 minute interval between each sample, __but what is the *actual* seasonality__? That's what we're unsure of. Let's take a look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD8CAYAAACW/ATfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wcmMxLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvjNbHMQAAIABJREFUeJztnXmYFdWZ/7/vvb3RQDc0NDvYgII7UZGIxH3DZWJmYhYTM0adYZxkkoyT0R9q1KwToyYmmTgaEo3RJCZq0BhRFhHEBUEW2ZodgW5k6QYamm56uX3P749bdW/durXXqVt1+76f5+mn7606VfXeOue85633vOctEkKAYRiG6d3EwhaAYRiGCR5W9gzDMEUAK3uGYZgigJU9wzBMEcDKnmEYpghgZc8wDFMEsLJnGIYpAljZMwzDFAGs7BmGYYqAknxebPDgwaKuri6fl4wUXT1JbN7XigF9SjG6pjJscRiGKRBWrlzZLISo9XOOvCr7uro6rFixIp+XjBSb9h3F9J+/jYlD+2PeHReGLQ7DMAUCEe3yew5244SAAOcjYhgmv7CyZxiGKQJY2TMMwxQBrOzzCIHCFoFhmCKFlX0I8CsEGIbJN6zsGYZhigBW9gzDMEWArbInoqeI6AARrdds+wQRvU9EHxLRCiKaEqyYvQNilz3DMCHhxLJ/GsB03baHAHxPCPEJAPcr3xmGYZiIYqvshRBLABzSbwZQpXyuBvCxZLmKgp3NbWjrTIQtBsMwRYDXdAn/CWAeET2C1IBxvllBIpoBYAYAjBkzxuPlehdqMM7FjyzG2WMGYPbXpoUqD8MwvR+vE7T/DuAOIcRoAHcAeNKsoBBilhBishBicm2trzw+vZJVu1vCFoFhmCLAq7K/GcBs5fMLAHiClmEYJsJ4VfYfA7hI+XwpgK1yxOndcDAOwzBhYeuzJ6LnAFwMYDARNQJ4AMC/AvgFEZUA6IDik2ecIYRAMsnLaBmGyR+2yl4IcaPJrnMky1JUJDlnAsMweYRX0IYEG/YMw+QTVvZ5RLuCli17hmHyCSv7kGBdzzBMPmFlHxJs2TMMk09Y2YeAANDDyp5hmDzCyj4kDhztzP7e2oFETzIkaRiG6e2wsg8BIYDLf/ZW+vsHOw9hyo8W4gev1ocoFcMwvRlW9nklFY6TSGZb8Ov3HAEAvLHxQN4lYhimOGBlzzAMUwSwsg8BnptlGCbfsLKPAOpaK8GjAMMwAcHKPgT0Op2UpbWs6hmGCQpW9iHAFjzDMPmGlX0IsKpnGCbfsLLPI2oiNE6VwDBMvrFV9kT0FBEdIKL1uu3fIKJNRLSBiB4KTsTeB+t6hmHyjRPL/mkA07UbiOgSANcDmCSEOA3AI/JF672wrmcYJt/YKnshxBIAh3Sb/x3Ag0KITqUML/20oL0rkfW9o6snJEkYhgmCZFKgo9tdvxZC5OiGIPHqs58A4AIiWkZEbxHRuTKF6k2saWjBqffPw/wN+9Lx9K2dxhXM7h2GKUx+MncTTr5vriuF/6flu3Hq/fOw+2B7gJJl8KrsSwDUADgPwJ0AnifSvocpAxHNIKIVRLSiqanJ4+UKlzWNLQCAJVvNf7t65wQ7eBimIHl+RQMAoM3EkDNi7vp9AIAdzccCkUmPV2XfCGC2SLEcQBLAYKOCQohZQojJQojJtbW1XuUsWAxHQIZheiUmNq9l2XyZeF6V/csALgEAIpoAoAxAsyyhehNOKpLdNwxT2Khd2I1xF0vnSZEsjAkldgWI6DkAFwMYTESNAB4A8BSAp5RwzC4ANwteFuoZvnUMU9ioXdiFYZ8eGPK17sZW2QshbjTZdZNkWXolTuqeVT3DFB9pN06eFACvoM0jZv48NuwZprDx8nSeZy8OK/t8YdUWWNczTO+AXHjt01F4ebL2WNkHjQsnHlv4DFOYiPT/6HZiVvZ5wqoJqCN7dJsJwzBOcGewFUboJeMQMvnMMEwvQmT9c0TGjSNdGkNY2ecRszpl9w3D9A7c+N8zxh/77AuKrkQSdTPn4Fdvbk1v+7ChBd95eb3FUSlUP19TayfqZs7Bos2pvHLj73kNtz79QTACM6Hy8LxNqJs5h9dYFBg/nb8ZdTPnIJnMrjezWuxJCtTNnIO6mXPw5qb96e2f+P58zK9PfWfLvsBQs9fNWrIjvW3xZmfJQPWVPWftXgCphvLmJk4o2ht5bNH2sEVgPPB/i1P1pl8IZTbvlkgm05//vmZv+nNLe3f6c5KVfWGijaXXh2GZWXFs2xUvbNgXFmarXdPROLrdTuo3XxE8rOwlYVSp2kq0jLP30EAYhsk/dn3Ti+JmN06B4iY3BsPwuF6Y6OstrbA9Wfb5gZW9JIwqLNuNY16lUV6IwTCMc/Q9OV9JzpzAyl4yVoa909BLVv7FA0fjFCZmfVa/3Ymy53QJBYZRhbFLh2F6J2YGmX57viJtnMDKXjJu3lSjwtZd8cI13zsw7cJOfPZRmaAloqeI6IDyohL9vm8TkSAiw1cSFhN+3kjFup5hCotcN47xdidunHz59Z1Y9k8DmK7fSESjAVwJYLdkmQoaL7lwWNcXLzzQ9xJMcuM489nLF8cIW2UvhFgC4JDBrkcB3AXWVQDM4uyt95vu4zvKMAWJ3iXrxGcf6dBLIroewB4hxBrJ8hQUB1o7UDdzDs747jx85rF3AQAH27pw46z3c8r++YMGmFXro29sCVJMJgLc+cIaXP6zt3K2c+RVtHl+RQPqZs5BR3dP1vZX1nyMuplzcKwzlSbFLBpHr/w7untQN3OOZZmgsH0HrR4iqgRwD1IuHCflZwCYAQBjxoxxe7lIs6bhCACgtSOB1o5EevvSHQcBcEpjJsMLKxvDFoHxwKMLUobYobYujBjQJ739sUXbAAANh9pxyvAq0+P1avzA0U7bMkHhxbIfD2AsgDVEtBPAKACriGiYUWEhxCwhxGQhxOTa2lrvkjJML4R99oWBvppiStRdbkK07HI5+0N8knNt2Qsh1gEYon5XFP5kIUSzRLkKArvHL46zZ5jCxqwLp18WLnT/beLsDVVGVCZoieg5AEsBTCSiRiK6LXixCgO3deTUimNjj2Gihd6wUw051XI367P6vPdG0Tn5svZtLXshxI02++ukSVNg8CM4wxQnlHbjpL6n89nbZsXMhfPZ9wK8rKZligs2GAobvcVvF2dvGKLNyr4QsK4lu4bg9DiGYcLBzGCLpd04qf+ZFbRe4uyjs4KWMYF1MuMXjrMvDPR9XU1f7tayNzL52LIvAOzqiN04DNM7Ib1lr/uv4sTSj3KcPaPgdkTmJwFGD7eJwoRM4uz1qtvRK0cjlAiNYRiG0aCPs4fJ95w4+xDddqzsPZBMClz08CK8smZPIOfftK8Vlz6yOJBzM8Fz69Mf4GfzNzsqKwBc+shivLSa0ylEETNPbP3eowCAG3/zPj5saMna92FDC+pmzsGk783HVT9fkrXPyIi/728bcMYD8/CLN7ZKkdkMVvYe6EwksetgO+Zt2O/qOKNRfWhVec62TftasaO5zbN8TLi8uekAfvnmNkdle3oEdjS34dvPF3VOwYLm8cWZuhaa70eOd+eUNfPYtHYmAk+IyMreA07nXfXlegxmZ9668xIJEjGFSncyCQAoiXFXLFS0ClyITKSOEWG+gJxbmAecKntHkzNMUaMaAPEYR25FGau+q7XhBASiOm5HVKxo41VpGx3H0ZnFTXePatlzQ4gizvpnpmPbWfZhGnys7AMkx41jUNMx1vZFTaIn1SZirOwLlqTejWNRlezGKRKMKpq7eHGTSLJlX+hoF04JCMvFlEYGX75gZe8B724ctuyZbLp72GdfCFjFx+vfNW1VlWHmvWJl7wGnCyP0vjujpdKs64ubBCv7SJPJgWNeRr/PyoDLVzpjI5y8vOQpIjpAROs12x4mok1EtJaIXiKiAcGK2TvQv8gA4Pw5xY4aesnKvnDR92qrmjTSAfnCiWX/NIDpum0LAJwuhDgTwBYAd0uWK9I4fRLT6/EwR3UmmqiWPfvsC5csn72ApbaPtGUvhFgC4JBu23whREL5+j5SLx1ndNhlwGOYRA9b9lHGyYN31qIqCJvQywJ64bgBtwL4i4TzFAxOqus3S3bgJ3M3ZW0La1R/dMEWlJXE8PVLTgxHgF7K0+9+hD0tx3HK8Cos3tyEV9Z8jPG1fdP762bOwTknDMTVpw8zPUe30ii2N7Who7sHFaVxKbLVf3wU331lA565bYq0cxYzVl33nW3N6c8vr/4Yf11lnOdICGEbjbPvSAeGVVd4EdEWX8qeiO4FkADwR4syMwDMAIAxY8b4uVxkcDI6/+i1jTnbwoqx/cXCVIIlVvZy+e7f63O2bW/Kzmm0ctdhrNx12PQcWh/uqt2Hcf74wVJk+97fN2D5zkNYvbsFU8cPknJOxp6n3v3IdF9S2Bt8zyzdibumnyxXKAXP0ThE9FUA1wH4srDQfkKIWUKIyUKIybW1tV4v1ysIc0EFE020kV3KXK2c8yqn5fn/6CBgrwOCjMP3ZNkT0XQAdwG4SAjRLlek6OO1OljXM3q0Cl5mR1cHEV7H4Y9M3nr/dSOEsD1PkDrCSejlcwCWAphIRI1EdBuAXwHoD2ABEX1IRE8EJ2L08FohbNkzerQtQmb7UBOs8byvP2SGRgvYP70FGZppa9kLIW402PxkALL0ejj0ktGjVfAyO3oy7cZhbR8ZhP3TW5A6glfQesFjhRjls2eKG23fT0hsH+qZ2LKPDgL2bpwgn/5Z2XvA63skOc6e0SMCsuzV87JlLwcZNSMcROOwsu8lsKpn9GjbhMwJWlVpsGXvD5m3L6Xs2bIvKHiClpFFls9eYvNQJwI5Gic6JIVwYNkHd31W9nmEffaMHm2TkOrGUf6zrpeDDDstFY1jF3rJln2k4Dh7Rhbazi3TGEj77Pn1OJKQEWfvwI0jcWGdHlb2Lnl88XZc9tPFno69/2/r7QsZsKflOO6evTb9vlKvPP9Bg6/ji53fv7cT8zbsA5DKYeKFxsPtuHv22vR3bd9f29jiSz4taZ8993B/KGPlnpYO/PcLa3yeTGDO2r2WJYJcQctNwSU/mbsJh9u7PR17tCNhX8iAu15cg+eWN+D9HQc9HZ8+z1/X2hdiTHnglQ34t2dXAgDe3HTA0zn+31/X4rnlmUFXa+l19ciPs2efvRzumb0OL640TnDmFCGAhTbthidoi5x0nhN+JI8MPR6ft/V1KILy2XM0jlRkjJlOajfUdAlM+HBSq+jh1b+ur8OgQi8zp+JG4wf17sl4QnIy+cqWfZGjLuLibhsdvK521SuNoNIlcCyAXGQ8ITmpEw69LHLS+oC1fWTwaoHlGIia0wSxqIrVvhzkWPb2ZdiyL3IyeU5Y20cFr4FRVpa9zNBLXsAnFxldz0mdBJn1kpW9C8J6M3wmZpqJCl4naPXugKBSHKvisc73h5pbKF+GFlv2EUFmVkI3ZCZoWd1HBa+Wvb4Og7LsGbnIeCG8Ez3ucymNJU5eXvIUER0govWabTVEtICItir/BwYnYnRIBLm8zQJe+h49ZFn2WekSZObGUTQLDx9ykGFoOcmW67VdOcGJZf80gOm6bTMBLBRCnARgofK91xOeZc9unKjhdTI1N86effZRJhN66f9cTqokSB1jq+yFEEsAHNJtvh7A75XPvwfwGclyRZKExBWObmDLPnp4nqDV9bigXl6inop1vhzyFY3jNyWKFV599kOFEGqSh30AhkqSJ9KE5sYxWCBzoLUDD8/blDNpvH7PEfzzU8vx2jrrHBxGvLZuLxZ5TAPQGzjS3o0fv74RCZsON2/DPsxXcuS4Re8OeOCVDenPdtcFgOZjnfjJ3E3pp4AVOw/hz8t355RT28zizQcwZ+1evLS6Ee9ta/Ykc7Gy7cAxbD1wDACwbs8R3+dbsUtvM+fSHaBBafsOWjuEEIKITCUkohkAZgDAmDFj/F4uVMKaQDOy7O96cS0Wb27Cp06sxdTxg9Lbr/vfdwAAS7Y0YeeD17q6ztf+uAoAXB/XW/jhnHq8sLIRp4+oxj9MGmFaTs2P4wUr+9CJZX/P7HWYX78f540bhIsm1OKGJ5YCAL44Rd+3Uuf68eubsrYWa9164XNPvCf1fGeMrMbqhhZLC39gZZnUa2rxatnvJ6LhAKD8NzUHhRCzhBCThRCTa2trPV4uGoQWLWHgs1cf9ziCQx6didQ9Dcvf7cSydyoju2/8097VI+1cU8cNQkkshvM1hpme8pJYoLmMvCr7VwDcrHy+GcDf5IgTbULy4mgs+0xLUCf6eCKu98ADd7SQOUdWEid0J5OWg3A8RoFGTzkJvXwOwFIAE4mokYhuA/AggCuIaCuAy5XvvZ4gc01bkcl6mUFtiKwe5BH2vXQ1QWtTNOzf0huQmWU2HiP0JIWlso8RBfpEZuuzF0LcaLLrMsmyRJ6wLC8j61218tmy7z04ifbiiKz8IdWyj8XQ3SMsY+1TLhxeQRsJwlKs6mW1V0/79ljX9xpkRnsF+S7TYkHmuFoSI/TYuHFisWAte1b2Lgg7GkfbgWNs2fc6wlq0xxgjMz1JSZxs6zdOIfvsmQyhKXuDpe9qM2T9IA8ra1iWpWx1FjeL9pwsvWeiQ0mMkOixrjUiCvSJjJW9C8K2orWXZ599cBhZdPkY6MNatMcYI9WNE4+l2pBlNE6wXllW9i4Iz7JX/2euH4/lbgNSFgQjH1lVb1U7btqXXaQImwASkDpBS+juSVo+kcUDjsZhZe+C0NZUwciNo1r22WVjrOw9Y1W9+XiCcuOzZzdO8Mi17B2EXsaI89mHyYL6/VjT0AIgPJfJlv2p/Bw/mrMRjy3ahkRPMp1QSxXpQGsHnlm6E3GNC6IrkUT9x0fzLG3hUf/xUdz70jos/yiVu+TtLU05ZWRV/bKPzPOjtLR3o6W9CwBwqK0LT73zUc6Tm1MFZCbv7FWN2NF0zOFZiou1jS1ZOY+kTtDGYjjY1oUVuw6bllHj7H+zZAdaO7rxYUMLFm7cL08GaWfqpfzrMysApHKKBOXGuWv6RCyo348R1X0wxyKB2YcNLfiwoQWD+5XlrKC9/dmVWLW7Jav8pn1H8elfvRuIzL2Ja375dtb3F1Y24uHPTcraJmugb2rttNx/9+x1ePymc3DnC2uwcNMBTK4biDNHDZBybQD4r+fXoCwew5YfXS3tnL0Fta/Izh8048JxWbpj7OC++Ki5LadcPEZ4Z1sz3tnWjK0HWvH8ikap8rBl7wInryW86bzcZG+r77sCOx+8FtNOTOXFePa2KVn7v3bxiXjpa9Pw2JfPdiTH8a6e9IIPVQkdbOvKKRdkBr1iw8md/Or5dXjm1in2BS042tGd9V/NheMWq6iOriBfh9SLkGXY33nVxKwnslED+9her01iXh4VVvYucJIuIW7QQtRNqjXu10gkIssXV6twrhV5OAmJi8dIyuvrAE20lUkdchBW8Mhy4sSIHM2laXWHkR7xLYf0M/ZinChPIz9f5h2y8mTRW/ZGUXscypdfZCh71SBQT6M3MPg9xPlD1r2OkT6vlfF5tQZcEHEWrOxd4MRv66R9yDDK1IahimQkW1hv1uqNOLmTUpS9cri+fu3g9AjRhYh0GWvNymU+BxFVx8reBU5cnUavL9N3QxkdM2PZq//ZjRMkTqosTvLcOOp5nNahXj6u+WjhxAjUth1244SMk45n1NeDsLr0uXGMBqIg32dZdDiowhjJ76ROo4B4JbV8ZNakUXpyPVnKni37cIlSh1LbQjpvDlv2oROPxSS4cVSfvTs3Tk5Vc9X7Rua4HXPgxsny2bOyDxevylP/pimZPntVJKNzchZFeThZsRqPybPI1PM4NTD08nHNy0CewtUOHOYTtJnPkXPjENEdRLSBiNYT0XNEVCFLsCjibILWvJLSe3yHXuZG4xgNRByNIw9HPvtYzHduIvUJLR2N49Fnz/hHbvScvWUfWTcOEY0E8E0Ak4UQpwOIA/iiLMGiiKPQS6t9ATQejsaJDjIte/2Tmx59devrn6NzooUTn712QDAK9PCLXzdOCYA+RFQCoBLAx15P9M7WZsMlxGHw5qb92NNyPGvbO1ubnSVCM5yg1X2X8JCtXual1Xuw62AbWjsSOWXe35Gbh6UnKfD8Bw1I6CZvP9hpnrOlWNl9sB1/Xr4b3/97Pf78QYNt+ZikaJym1k7Mr0/lRPn5G1uyFLf27Frj47nlDXhzU+qY5R8dCmQFZiGxaPMBNB5u93Tsgvr9WLbjoG1qCzdolXdLe7dhmXiWss9s37yvVYoMnnPjCCH2ENEjAHYDOA5gvhBivr4cEc0AMAMAxozJTSWgctOTywDIz0vhhVufXoGBlaVYff+V6W03PbkMD332TNtjayrLcrZV9Und5vSbBAXwidEDMGKAd6+X2nhW7jqMix5ebFjmr6sac7Y9v6IBd89eh8PtXfi3i8ant3/uiaWeZemtXPjwIlflzzlhIIZW+fNkCgE88db29PdN+1rx3vaDmHbi4Jyyf9EMQD94tR5Aqv98/tdcl7f87gP0KY1j4w+muz5WzYfllfKSGP7jkhPx0wVb0tu0hvqKXYdxy7Q6zF2/D3uPdKS3m7lxrvr5El/yqPhx4wwEcD2AsQBGAOhLRDfpywkhZgkhJgshJtfW1nqXNM8cNhh9jdIlXHPGMPSvSCnzVfddgYrSeHrfg/90BnY+eC3KS1LbtI9pL399Gv7vy+d4ko3gfYXdISWHjtHvY7xz8rD+OGvMQFSUxvHuzEs9n0dA4OCxbIuytSO3rgSAw+25+ZCYDMe7g3+62fngtTlPc5t/eDW+cdlJWdv0/fWBfzgNS+++LLuMplDUonEuB/CREKJJCNENYDaA8+WIFR5Wvk4jnz2B0hOuMcoewc38bjLcqX6XcvOqe7mUxOUtddfXrbbZcb1FDyf17ba/Ri0aZzeA84ioklK/5DIAG+WIFR5W4YqG0TikCa3Uvz+IcopKw2tbUAcz1hlykZXEyqiJadsdz7tGD7u3hgHu+2ukLHshxDIALwJYBWCdcq5ZkuQKDasIFmPLXqNAY8iq1SAVqleFItJPIazuZaLtnH6euoyVfe42rr0I4cSyd1ljQVj2vl5eIoR4AMADkmSJBFax6UbKXr8yTjsgm3V6GcaZ17ag/gTW9XLRxtf7McoERI5aMEpzzAZ+dHBS3U7ahNE7pmXCK2h1qJZ9jHL990ZuHNK4cWKU7cjR16+qYH3HQBvks3eKGvbJul4u8Sxl79dpn/1V2+54kI4eTuqkoN04vRXVZ18Si+U8PpvlFdPmq89eFq0vKS9dgldXQVpvsNaQSkks05X8dFQhch/5OetFtHHionFrAERxUVWvQ3XjxA3e9G5o2UNrLWdXu76+gnh5iVsyTyHSRGGg99l7P4+RXneTgI9XzuafIPpSEN2Tlb0O1Y1jpOwNJ2iVN8KnPut9+N5DL+0akOcGJjIDEyMPrc/e753VDxZGCtxMqbOuzz9OnrKdlNHWXRBPc6zsdahunHiMcl71Z6zsddaYhRsn89VJXnzrxuFVWbMXJxjikqJxIHIHi2z3ofW5WdfnH1kTtFpkpFTR4ysapzfSo2j4I8e78eQ7O7L2LVRyj2jRLqpKfdd89tHpU8reuMKX7TiINY0tns6rDljr9hzBe9ub8fbWZgzTLfF/cWUjzhhZjYnD+nu6Rm/gx6+5WzISJzmWvYDAJl0ulKQQeH5FA06oqUxvm71qD04a2i/3eA+m/YHWDry9pRlnjqrGSUN7d50fbuvC1gPHMGVsTXrbos0HcNqIKs/ndDRB6/KcW/cf8ySLFazsdWgXVT0yf0vWvvV7juaUv27ScIwc2Ae/XLgVpfEYSjUxU5NGVbu+/uiaPuhbVoKrThuGXyzcaljm1bV7XZ9XZcv+lCJZUL8fC+pzBy8A+O8X1gCIRp6isPj1kh32hTRoLfvyEu8PzEKkBmItTa2d6bZw+SlDAQBzN+zD3A0Gx3u45tU/fxsHlTQavaHOrQa8L/12GTbuPYqPfnwNiAjrGo/glt994Ot6Rw2SEALAJ8fWYOPelM5Y23jEsIzK4H7ZObVeWr3Hl0xGsLLX4SQF/B2XT8Cjb2zBNy89EZdMHIJLJg7Bf10xAQBwhqLg4zHCCYP6Zh2XCb00P/fbd2XyqtxxxQQc7+pBSZywctdhfHHW+45+w6RR1bjq9GF4aO7mnH35yBdSlGhMt5J4DDsfvBaH27pw1g8WAEjlUHpt3T5Pp27rNFYmRnjx2auKvrdgdQ9U5StEqj96yS+088FrUTdzTvr7RRNq8daWppxyf/m3qenPR45b56Ja8Z0r8OXfOuvfXmFlr8ONr8yopNrnjawL8hB62acslUTNlbVIOYkb0nQn2KsbBEZ3W/t471QJGxUzSsBnRpRenRkWTu6BzLkrJy+s8ZvLSgY8QavDSV+xstDT7w61OM4LbhqLfiWvlm5+e1Xe0NaZY2VvUNBoBS1jjpO7JXNQdLKuwkn3DXqcZmWvw5Gyt9oXUKW6GSf0i7u08NurgsFoMM6y7B0+z/m17P0MKr0FJz9NprJ3kscmCutaWNnrcNIp05a9Qdm0ZW+UvMriOKfXdFQW5qGZXQm27IPA0I2j+ezHONeGXtq1A8eDSu/V9Zb3wMm8mVucvJ0sCutaIqHso2RlOOmU+ve/OsVPhbs5NkZkqhTYjRMMRvc75smNk7utx0Wd+Zkb6C04uQcyVY4TN07MgaYtCjdOlFySbgYeo5JBJDACXFr2FmW7zRL8MNLJrgfvbpxujevNrnk6bb1RMrDyiVolct04Tq7Llj2AaDU8NwOPybtMPB1nhzs3DplO6HI0TjAY3W0vlr1Rwa6sAdr6RE77Um9uBVa3QO0X+Z6gjYCu96fsiWgAEb1IRJuIaCMRTbU/KpcoWfZOuoETn73xgV5lcmkZkPmlrPL1M96xi5by08S7NfMsdn3Fx5jSa7D02Sv/ZeocJxO0UTBo/cbZ/wLAXCHEDURUBqDS7gAjgsgD4RVHPnsLxesoGseFPG7Omy5rUZ4naIPBzrJ3ilHb0K7qNsrPlHV8dLpSaDi6B3meoLWrNyB4PehZ2RNRNYALAXwVAIQQXQBe/pvWAAAgAElEQVQ8LcULu4HuP9qBrkQSo2sqsXLXYdvy6T5sEXFjeJx6mIcf7NZnb27ZO7/2wWOdGNSv3PmFI0ZXIomNe49iaFUFEskkRg30ZIt4JntRlbP7brSs/p1tzenPdu6Hd7Y2W+5XWbnrMI4c74K+pQghfC8AOtzWhUPtXRhfm5u7RwZCCKzafRhnjxmYI+uKnYeweX+ryZGZOvlg5yFcdsoQKfI4ceNEYarMjxtnLIAmAL8jotVE9Fsi6qsvREQziGgFEa1oaspdUgyEr+w/+T8LccFDi7Bw4348+Pom2/KTTxgIAJg6flDOPmurPz/RODub28199i5a3QUPLXJcNor8z2sbcf1j7+K8Hy/Ep34S8G+xicbx4zbQPo3Z9ZWv/2mVo3Pe+Jv3cfsfVuH2P6zM2r69yX8CrisefQuX/fQt3+cx44UVjfjs40sxd312+onuniRueGIp7n1pve05/uWZFXhlzcdS5Ll4Qq1tmStPHZr+PHVcRm/841kj058Pt1mnVPCLH2VfAuBsAI8LIc4C0AZgpr6QEGKWEGKyEGJyba3xTYnKEu9dB9tty9T2L8fkuhrUf/8qXDwx1zKwGuT92EtuxomjHd2m5Z08Tqq0dxV2Hh19QrEguO+6UwEYD8baLVb3fcsPr0b9969ydL2g+0pbp/86bz4WbK4dNZlfw+Hs/uqkbWvradfBdimG5pWnDbOtv89NHoVV912BD++/As/eNiW9/ZHPTcKmH0wHAFw00X7Q8IMfZd8IoFEIsUz5/iJSyt810VD1zuSoKE3dssoyYw+Y25cUOMXNQBGPmUfjRGsyPFjyEQDRv8LcE6qtAquJ8bKSWNZrDa0IWtm7cfOFhbqiWD8nEqbRaKYPVIgINX3LMKCyDCWazLjxGKGiNJX/qtRJDKcPPCt7IcQ+AA1ENFHZdBmAei/niopl78SvaudO8ZtKQcax8VgUonqLBKXJGNWPdsC1C4JyWr9B62I3T35hoXZT/cSoIzUS4Y7hJKrHD36jcb4B4I9KJM4OALd4OUlEdL2rJGhmBPGiYOXKjkuWxMxX0DJyybx/2Bq7kFfH1RVwX0lEYSbRhh7N2+S0ODEao9wtglqQqeJL2QshPgQw2a8QUYhBBZw1Fltl7ij0MthonJJYLBIr9ooJu/qRZTGzG8fcjePIsI9wtwjaso/ICtqwJUgh4ynQyQStl9/rZkCMs2UPID8d22m12ClRp5FawSv76Fv2SRPLXjgQPcpGUNCWfTSUfdgCKDh6DLQz7C0K+Am9dPN0XRKPcpPuXagtxu6O21n2Tusr6AzVhZACO+3GybHsoy+7FU4WZ/khEso+OhO09mXsFLaT+vLyc924AdiyT5GPIU+tS7v7bW/ZO71esH2lECZoVTeO/p45y1gbgECSYDdOHnHSkeyUueWiKvU6LmRScTMglsTMX0vIBENv8dl3F4CyN3XjREWReKQ43DgRqSRnb6myCb20jL1Ur+P+97qZOIvHYtEOO8gX+fDZOxy6baNxnPrsA3apu8mdb0dQr1NUPU250Tj2x0a5WwQcZh8NZX+8O9iVmomeJN7fcRAdBtc53JZZ7fe6bvm1F4J6EnNjGRKCb9SH2rrQ0h7sSslCIDNuW99xWTp0R7P/dAZWtLSnlux/3HI83V8aDrVbptk43NaV1Y9UgnoKMbLsGw61473tuXmBGg+3o+FQO9o6E9h3pCNr38ctxyOVGDBon73fOHsp3Pb7FYGe/+F5m/HrJTswYWg/zL/joqx9Z/1gQfpz/d6jtucKK85+aJXzhGREQcb7pzhbuW87H7w20Ov4IR9WXHqC1uJi542rwYSh/fHM0l0AgIsn1mLxZuM8UXZ0dAernD7YeQi3TBuL8x98ExdNqMXPPj8JFzy0CF857wT84DOnGx5zlklb6BEiEAVjNEFrlsdJzYk0uqYPGg4dR7/yjER//qABf/vQXX4cbV4bPeMG98WO5jZX59NyyvAqAMA5Jwx0lJDRLZFQ9tsOBGutvLf9IABgy/5grwPYrKBV9nqxd0YNrMSzt03BV55c7kwOnSBL7rwEFz5c2InNgmJ8bV9sb3LWSW+dNha3TKtDZyKJIVXlaWVhVu+r7rsClWVxlMQI/3rBOPSvKEGfsjiOdSRwzg/fML3OzVNPwO+VwcEtZ40ZgNW7Wzwdq132/9aWJhw5nrL0tZk3nRKUyykzQet8OG84dDx1jG671qtw45TReG55g+HxS+++FGXxGPqWm6vMOd+8AJ0J716KyXU1WHr3pRhWVYH6vUdx/a/elbruIRJunKCRGZJlZzFb7fdrbDtNGat/B21N3zKMGWSf3ndwvzKvohU0eqUxqK/5fajtX47RNZU4cUg/VFWU2k701PQtQ0VpHCXxGEbXVGJAZRnKS+K2qaOr+5Q6/wE6poyt8XTcmJpKJHqShnNKXuaZgnbjeOpPFsdU9zGv9+HVfTCoX3k6j40RfcriGFDprw8Nr+4DIsJpI6pxxqhqX+fSUxzKXmKbs4+zd3ASj/I4bdypfPbue0JE5snzjl4pWflO9YaDEzeOJ3ycsDzurVuXxAmJpMia6PS1NiSgBqWeV3ZgR9TCMmWHwRaFspd5z+wse8tFVcp/r08aThU4wduLM3qbrnceu5793ctEmexQVz9zdWUlHpV9jNCTFFmDn59wYScrWr2gKkHZ2WMDnh91DSt7D8i0APKVudLX8VEzUSKOvlO5UfZBPQ35GTy8Kvt4LIbuHpF1P9LvW/ayEDAoN45yXtmRnVFbm8LK3gNS3TgSzuFVHse6HtlPGL3NYneK086rd+OUuFL2PvzHFvixMstLzP3KVpTGCT3JpGH79PI0GtRqXHXi11tCQfMby5Z9L0BqzgwJrxb0LI0rnz3jFL1ys1rJqC+byY0jFz+Dh3fLXvXZa904SpsNOHmfG3oCsuyj9kTMyt4Dcn328s4VFLk++9BECRXnLwTxbtlnriW3Yfg5X6nXCdoYIdGjU/ZRdOOkffZyU4VHrW/LTjftW9kTUVx54firMgQKAqk+ez/H+vXZO52gpezcOFFJRxFV9MrezYK0wHz2vtw4XpV9TJmg9S8DEJwbJxON4/5YJ+tgokIULftvAdgo4TyBITf00n+D8Oyzd+rGQbaVUqyq3un90mcCKHGRpCSoe+tnBbTnaJw4oTuZTFvOMSJffSeogTBt2Uu++1Gz7COl7IloFIBrAfzWSfkeIbJybHQlkjjWmfAjgiOMbtmR492eXsEmw7L3Hnrp/DpedIXRE0BTayd2HWxL70smBY4o+VMApPOOqLR3JdDU2okj7d1Z5bR0JZJo7TDe5xUhRDpXz7HOBDq6e9KrP1V6kqky+jwu+t8dt3j5t75sUBO0fk7nJ/Ry24FjaD7WmZYhKcxdJkeOd2cppONd2atHg7DsO7p70jrjo6Y2bG865qof69uEloi57KW7cfymS/g5gLsA9HdSuP7jo/jGn1bjia+cAwD4wqylnpd1u0H/mN7R3YNJ35vv6VxOG8TpI6tytp06ohpAA06o6evp2uW61Xtqvg89JDEV2rk/Si3pf+izZ+Lz547GI/M34/8Wb0/vV3OSqHlR/vnJ5Vix6zCGVpWjJBbDuzMvzTnnV55chmUfHZKaV+fJdz7CD+dsxNt3XWKaJ+WR+ZvxuCL7Tz83CZ89ZxSA3PZx7gkDsaYh1S77V5SgtcPcIFFXNZ8yLLe+/eDUsu9fXoJWncHkdVFVWUkMrR0JXPHoEgCptm6mbtq7Epj0vfn46vl16W2n3D8Xz//b1PT3IHz2k3/4RlrZ//LNbfjlm9vwL58a6/h4K/25+1B7zraK0ljg+YjMmDp+EP6+xl3uHis8K3siug7AASHESiK62KLcDAAzAKBs2ImYuyGTWVKv6CcOdTRmuEbf5jodVt78Oy5EeUkMFz28OL3NiRvnrTsvNlwOf9Mnx+CcMQNx6ghviqGfLi/Hb/55Mqb//G0AwN++Pg3XP/auIqNuUNL4YNV78bPPT0JHdxL3vLQuq9jyey8DAEz50cKsa729rRmfP3c0Xlu311LGFUoCp/1HO03LLPvokOU5vPDGxv0AUk8aZryukf3trU0aZQ989uxR+NZlJ6Ej0YNxg/viG5eehEWbD+CSk4dg9e7D+OOy3VhQvz/nnJecPASvffMCnDJcbtslApbfcxlAQHNrF6755dtZ+5ffcxm6kwJ9y+LoTCTxwN82pPtWiU7ZnzeuBiOq+2D26j2W17z76lMwb0PmNxIoY9nryrZ1pqz4V9dmK6O1jZk+HcRckZEnQK17v7S0d2NgZSkOK0+kL94+FeNr+6E94Ky8Zjx8w5koL4nhxZWNUs7nx40zDcCniWgngD8DuJSI/qAvJISYJYSYLISwfTF5UO+/9JqjY/TASpwwKNsKd2JvnTCob45iBlIDhVdFrzJucEaeytLMNSaNHqC5jrGcI6r7pD//09mj8KVPjskpM6R/BYb0r8C4Wm9PH2HhJKxVG1KpffteUghU9SnBmEGVmDC0P0riMVRXluIzZ41EdZ9SXDxxCCYMTVnwRk3p1BFVgUTjDKlK1YVRmxlSVYGRA/pgQGUZhlZVYPrpw9L79L5nAqUzKloxrLpCf2D69zrtQtr74MFL6oluSa9STAqBoVWZezBiQB8M7FuGkQP6WBwVHBWlcZw4xFk+LCd4VvZCiLuFEKOEEHUAvgjgTSHETX6ECerN9jInREP362mubyYLgQwXVUXl9Y9B4KRetClxtS/WSCaFfRqMPEdquL2adiDTrxNw2mb1IaeEjHXudJ5Je4Z8veJQlpHYkxRZ/cZLCK5sZEoQqTj7oF52rH2cdPNoabRsPuzwLCdXN7PsXSn7Ah0XrH6itj61iigpgn9xhFvciqMdyPTvMnWq7PX3wMpn74R8GReyLPsekd1v9O6wMJBpXErJZy+EWAxgsd/zBOXG0TaFnqRwbKUYWXsWgRp5QWt5mFr2Op99OpLGRZ/Qd1Rh4ruNCk46hbY+tZOHSSHCf2LT4dYtpNVL+nabetJzf80sn71jN07mc/6UvRy9IXTtIAoGgEzjMvyhS0NguTQ0jc6Nq8iorsO27LUymSkEfedOu3Fc/PagVj8GRcZnby63tvNq74UQzqNf8nVX3OoZrfyyDBLthL7+d5s9IYfixpFl2evkjYQbR6IIkVL2+fDZJ5LC02RTZpskoTyiHWzMRNHns1d/r52lpd2tf8hS70X4zd8YJ8v6sydosy17u36d93p3bdmT4WfA37oOM8verKtq+0yedL1kn33mu5vFdYVAtJR9QD57baNL9CQLeqLS6WOyka5w0/kKNcWCldTavqu14nqE/QRt+vx5ui2uJ2i1lr3ut3iVOXsFbfZJnPShQvPZJ4XIMpJKwvbZQm7OpUi8g1YlKJ+9tqEG9fQQBma3S99AVMvO3rLPnrQ02lfIdy/LjaNZGerEjZNvG89tuoSsaBxJyt4q9NLM4ieTAbUQ0Fv2EfDiSJUh/KFLQz4se/2beNwiO57aD2aDo/bxW4ubn50PqyyIpwerc2qVYDonusjdZ3n+PA13bptZPMuylyQDMr8312ef+m/VTgrtCVqvfqLQ1ws+9PJwW5fhqJ9ICiSTAkIItHUmcLyrx3MOlWQy88Yd7bW6e5KBvfU+H2QpLJPOFKPsBS1qMTeWlv7cR453ozPRYzpg9Ch1l3MerbskmZ0b6WhHAkIYH2cok9I21M/JpEi/IFvtmJZuHI0W7OpJor0rkf6dtgoyzx3f9QStpidL89lr3Dgd3T1I9CTRkxTptgBYu1CSSXdBASr6Pq/ty0GSTIrIzUkVvBvnrB8sMN037p7XMKR/OQ60ZpbbL/7vi1E32N2KznH3vAYAeO2bF2QlP1q0uQn3vbzepcQZJkhc0eYFbd2rb7qv6VumK0OoKM30/kmjUqtrzxxVbZiqoDRO6O4ROHNUZhXu6SOrsXhzU/r721ubMfE7c03lGq/cbz3j7nkNOx+8FkfauzHp+9n5iLT5iZzkyRl3z2s4b1wN/jxjarp+AWDm1SdnOqmFTth3pCP9eeWuwzj3h2/g5a9PA5BS/laMHphaRTlqYKWtnDJwG/U1WJOeQ/+UMsFHGhJ1MGzv6sGJ976ecz19+oL7/7Yh/fk7L6/DzoPtrnIgvbCiAXe+uBYv3j4Vk+tq8P6Og/jirPcBAM/cOsXz73DChKH9URonrGk8YltW7TNBE7k4e9loFT0AHGzrcq3sVVbsylZuizYdMC37xE1n4/Y/rAIAvHB7JqHTwm9fhMt++hYA4K7pJ3uSQxaqgXPnVRMxamAlXvvmBRhSlep4915zCn702kYQgKnjBuG7/3AqkgK4YXIqB8xvbp6Muev34cKTatPne/1bF2Bwv3I0Hm7HSRql8KsvnY3TH5gnTe6mYx32hRzw/o7cwer5DxowZpC9Em7RZTxs6+pB87FUBszTRlRbHnvDOaMwvLoPpp04yIW01iy9+9KsjKGvfuNTuO5/3wGQ28kf//LZ+Pc/rjI914Sh/fHoFybhlOFVWf77Oy6fgNsvHodnl+7KOeb+607F91+tBwD87pZzc/afPrLKdOxUM2NasfNgKk+R9snLjve2HwQA7DrYjsl1NXhrS8bgWLTZvO8aMesr52DGsytty117xnBcNLEWn540AgBw8cQhGDXQOkXCu//v0nQOnSCRGesfSWWvR6Zv1+xcE4b2w/TTh6e/n1tXk/6sZjYEvKePlYUq/3njUkpHmzdFOyCWxGP46rSxWcdWVZTi85NHZ21Tc6bU9s9O3NavvATjavtiR1ObPOHzgFuXhWrRD+pXZlmOiPCpkwZ7lsuI4dXZCuX0kZkBR68cPzFmAOz4x7NSg7r6BNOnNI5vXX6SafkpYzNtfECf0pz9ZfGYlL7XkxSOwxjTbrVY9ncvXHnaMPtCAMYMqszqF9o8Q2YMqarAkKoK23J+kRnrH6kJWjP8uOvs2oqqvAslSkf1XRo1Atm51aPmv7RCldVqPsbo93QoGQ2jsFpSi14cN9E5qqK0G/icKGAZdpabvqUWVX+v9vpBLWiMWNVnITP8s0CUvfcWZ3dsmbLO3MtEUhiov8dIOam/IApRBPlG/c1WisXotqjKPgqrJbXoZXVTpUaK0gi730xEUhZGuVP2qsGSkk07MRtUJJSft4IFjcyFXUWg7LO/6xVhqXIzCyU9gCqm0Uul05Z9PgWKCOpvdhu1ob7bIAoLaLToFZA+uZkVcSOr2OB4u99MkONC7XExkZl5JaLyXbv2IyCDLMrGkcwnzmi1cBP8hEraNRDVjeOmQYaJpWWv/IQIt10pGCkgMlAOuWVyb0yHEkIY9aXx7tw49nmCAHtFIsuy73bRgfWLtbTV2R2QsnczkOYbduNIPDat7AvEslfbu6HPXvkfdrK2oLHq827bSsayj9Y90yt3V8reQZ4gwH6AS6U4ljNB6xS1qFqP2voMyiCLWNVnUYQTtO4q2cmyfxXVZ18oS7tVOYvNss+uU6O6Iot92hLZZHz20eoKOT57F+LF05a9Nc7cOM6va4Ybn73QKXmvGWvdoH/ZS5SIhM+eiEYT0SIiqieiDUT0LWlS6XDb4LRtQt/5c332haXsrSJuVCssyhNOXsl+2Yi5G8cq5YbhBG1E3Tj+LHvVZ28TjWPrxpGk7F3km09b9soh2kODyp0V5e4i0wjxE2efAPBtIcQqIuoPYCURLRBC1EuSLY1bRZylGJz67AtE2etD07SkO2YEG69fpaG16ozO5WSC1si91RFRN45eGjd+5bSytzgfYD/AxYh8uVBV3Fjkav1pE9WlzxOYGydada8lEhO0Qoi9QohVyudWABsBjJQlmBa3DU5rAejfYqP/ng69LAxd7yj0Mmp0JZK+76+lZU+ZerVKe2BkGbZ3pVawRuEVdFpy3hrlos+n3Tg299xOkfQkhZS5LL2S1lv62rrtSqT2GU3QdiaCseyjPEFbKvGJU8oKWiKqA3AWgGVezzG4X7npEmz9kufLTxmC3958Lv6+5mN847nVWHXfFen8MG/U78e/PLMiXfaXb27LOvbtrc1Z3/uUpfLLnDLcOn9IVUUJjnYkLMvkg7rBfXGgtTM9SGlRV0GOHGC91NspYwf3xXYJK2gnfOd12zJ1M+fgrukT8dDczbhlWh3G1fZL5zD65qUnZtXjqfdnp3HY0dSWXumrzc2iR02NoOW55Q0A5HYqGehdMG4sPLXo+NrMimqjFcLac/Ytz1UF8+v3Y379fsfXNUM7yD67dCfuM6ijG6eMxowLx2PpjlS6hLv+uhZ3/XVtVpk3NvqXxYj+FdFNJBCpdAlE1A/AXwH8pxDiqMH+GQBmAEDZsBNNz/P960/D1yxyf2h5Y2MqR8bT7+0EAOxoOoaavqml30u2NpkdZsjIAX3w8A1n4pNjU+kH/v4fn8pKIqYy744L8VFz+KkDZn3lHKze3YKBfXM77wUnDcavvnQWrjzV2TJxO372hU9g8eYmLKjfj9IYYfbqPVLOa3q9+VsAAL97dycGa5STfsB2yrVnDMecdXtty02pq0H/itx0AWGgJgHUP6GUxmN47EtnozPRY5vHh4jwzK1T0qkwAODTk0Zg75EOnDmqGoP7laMrkUR5SRyPf/lstBzvTidLW/TfF+OSRxZL/U1ay33lrsOGZZ5b3oDrzhzh+JxfPb8u3f9VFtxxIer3HsX42n6GaU2mnTgI37n2VLy9tQmnDK/CV55cDgD49CecXzffjBzQB3deNRH/8RP/5/Kl7ImoFClF/0chxGyjMkKIWQBmAUD58JPStf7Dz5yO72iyT/oZXf08aMZihM9p8mKcMcq4Iw2v7pOTyyQMBlSW4ZKThxjuIyJXHcaOqopSfHrSiHSCKFXZnzq8CvV7c8Z132jrUcYcysyrT3ak7KPU2c8fPwgvf/ix4e+/9szhBkcYc+GE2qzvRITbLxqfU+7qM7LPOdYk4eCjX5iEO/6yJmd7PEboSQrU9i9HU6vxk7k2O6RVrbp5cbg+l9OkUdU4aWj/rGR+eq6fNBKnDK/KGgQBoLIsupb9kKoKfP2SE/EfEs7lJxqHADwJYKMQ4mduj9dPiBmtCM0HEZuXKwjy8VIKGZNxTkPqouTCUWWOWsBAeUnccLvaj43ciipOf4ub35ybTsJBHUanmkPBj4adBuArAC4log+Vv2ucHqyfEPOTTdJPHUZ5Jj6q5GP9mYyYaqcDeTxCMfbqZGHU3vJk5NoEMkZaucl+INtnb1UlbupcP6nKRps9np9fhBDvwIeezbHsQ+pwrOzdE5QiynoJuBRlX3iWvRoO6cKjkRcqTCz7uAPL3ulTmpunOX3Vcj+2JzSTRh/jG9aiFm4k7smH1ekmn4oZTus2SumNVZmjlr6jvNRY2asDpdU9dDpwO1k0pV7Gz6KzYiU8ZZ/js/deWb4maLmNuCYf7mQZus5p3UYpVYKqtKKWctvMjRN3MMegnXi1+llOLHv1/uQo++hUYWQJ7Rbp/aR+/KZ+JvOinBcjqkTNn2yGU2svSqtnnSjPMKgwsezVgdKqSTh1zzn5zRllb7ydMSc6bhwfHc5Px+A24p7epuzjEfLZxyI7QWui7OP26ZS1E69WrhonE7RqlerdRlFyxUWVyLhx/Pjs/fh32SJwT0D5qKTjNFNkWMEBRqjznFGz7MtNouWMXpSiR6vgrZ7CnfnsU9eL8gtHokqIbhx5I7M2z7VbgyjKeTGiiswXwAeJ07qNklWYjrOP2D02s+zVeDyrJxGtgrey3p24Y9W64gla94S2dEy/iKo8btKYTDjt/rlo60qlp/3an1ahVGkEbpMlsc/ePbX9y/HxkY6wxbDFqRIvK4lOG6ipTKWJ6GeQqyZMzCz7YVUV2NHUhmHVFaZ5lO59aT2++0oqH85x5f0BRjz4+iZbOQb3K8OxzgT6VZTgpCH9sPXAMQC5K2qNsAoPLQYon1YaEbUC2Jy3C3pnMIBm21LhUggyAoUhZyHICBSGnIUgI1B4cp4ghKi1K2xFvs2HzUKIyXm+pmuIaEXU5SwEGYHCkLMQZAQKQ85CkBEoTjmL+7mGYRimSGBlzzAMUwTkW9nPyvP1vFIIchaCjEBhyFkIMgKFIWchyAgUoZx5naBlGIZhwoHdOAzDMEVAXpQ9EU0nos1EtI2IZubjmhayjCaiRURUT0QbiOhbyvYaIlpARFuV/wOV7UREv1RkX0tEZ+dR1jgRrSaiV5XvY4lomSLLX4ioTNlernzfpuyvy6OMA4joRSLaREQbiWhqRO/lHUp9ryei54ioIuz7SURPEdEBIlqv2eb63hHRzUr5rUR0c57kfFip87VE9BIRDdDsu1uRczMRXaXZHqgeMJJTs+/bRCSIaLDyPZT7aSYjEX1DuZ8biOghzXZ591IIEegfgDiA7QDGASgDsAbAqUFf10Ke4QDOVj73B7AFwKkAHgIwU9k+E8BPlM/XAHgdqbWC5wFYlkdZ/wvAnwC8qnx/HsAXlc9PAPh35fPXADyhfP4igL/kUcbfA/gX5XMZgAFRu5cARgL4CEAfzX38atj3E8CFAM4GsF6zzdW9A1ADYIfyf6DyeWAe5LwSQIny+ScaOU9V+ng5gLFK34/nQw8YyalsHw1gHoBdAAaHeT9N7uUlAN4AUK58HxLEvcxHR5sKYJ7m+90A7g76ui7k+xuAK5Ba7DVc2TYcqTUBAPBrADdqyqfLBSzXKAALAVwK4FWlUTZrOlj6vioNearyuUQpR3mQsRopJUq67VG7lyMBNCgduES5n1dF4X4CqNN1fFf3DsCNAH6t2Z5VLig5dfv+Ean3UOf0b/Ve5ksPGMkJ4EUAkwDsREbZh3Y/Der8eQCXG5STei/z4cZRO5pKo7ItdJTH87MALAMwVAihvp16H4Chyuew5P85gLsAqPkfBgFoEUIkDORIy6jsP6KUD5qxAJoA/E5xN2ktOLAAAAMISURBVP2WiPoiYvdSCLEHwCMAdgPYi9T9WYno3U/A/b2LQv+6FSkrGRbyhCInEV0PYI8QQv+29CjJOQHABYrL8C0iOjcIGYt2gpaI+gH4K4D/FEIc1e4TqeEytDAlIroOwAEhxMqwZHBICVKPpI8LIc4C0IaU6yFN2PcSABS/9/VIDU4jAPQFMD1MmZwQhXtnBxHdCyAB4I9hy6KHiCoB3APg/rBlsaEEqafO8wDcCeB5IvmZ3fKh7Pcg5TNTGaVsCw0iKkVK0f9RCDFb2byfiIYr+4cDOKBsD0P+aQA+TUQ7AfwZKVfOLwAMICI1xYVWjrSMyv5qAAcDlhFIWRSNQohlyvcXkVL+UbqXAHA5gI+EEE1CiG4As5G6x1G7n4D7exda/yKirwK4DsCXlYEJFvKEIed4pAb4NUpfGgVgFRENi5icjQBmixTLkXqaHyxbxnwo+w8AnKREPpQhNeH1Sh6ua4gyYj4JYKMQ4meaXa8AUGfeb0bKl69u/2dl9v48AEc0j9mBIIS4WwgxSghRh9T9elMI8WUAiwDcYCKjKvsNSvnALUIhxD4ADUQ0Udl0GYB6ROheKuwGcB4RVSr1r8oZqftpcG0n924egCuJaKDyBHOlsi1QiGg6Um7GTwsh2nXyf5FSEU1jAZwEYDlC0ANCiHVCiCFCiDqlLzUiFZyxD9G6ny8jNUkLIpqA1KRrM2TfS9kTJCYTEtcgFfWyHcC9+bimhSyfQurReC2AD5W/a5DyyS4EsBWpmfEapTwBeEyRfR2AyXmW92JkonHGKZW9DcALyMzeVyjftyn7x+VRvk8AWKHcz5eRimCI3L0E8D0AmwCsB/AsUhEOod5PAM8hNYfQjZQius3LvUPKZ75N+bslT3JuQ8pvrPahJzTl71Xk3Azgas32QPWAkZy6/TuRmaAN5X6a3MsyAH9Q2uYqAJcGcS95BS3DMEwRULQTtAzDMMUEK3uGYZgigJU9wzBMEcDKnmEYpghgZc8wDFMEsLJnGIYpAljZMwzDFAGs7BmGYYqA/w8P5tUhXO5cnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# extract the data we're interested in\n",
    "n_samples = data.shape[0]\n",
    "xlab, y = data.time, data.occupancy\n",
    "\n",
    "plt.plot(np.arange(n_samples), y)\n",
    "plt.axis([0, n_samples, y.min(), y.max()])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__We can see that even though the data is separated at a 2-minute interval, the seasonality is VERY DIFFERENT (and much less frequent)__. Therefore, part of the issue is the assumption that the actual data frequency is equivalent to the seasonality. *It is not!!* It looks like we have 3 seasons in a 2 month period, which means we likely have 18 seasons over the course of the year (3 * 6). Only 3 are in the sample, however, and the algo is time invariant (i.e., it has no idea whether this is a year or a month; it only knows there are `m` seasons). Therefore, we'll set `m` to 3.\n",
    "\n",
    "This is really important to understand: the `m` parameter must be apriori knowledge!! Even R's `auto.arima` will have a tough time fitting a model for a `ts` with unknown frequency."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stationarity\n",
    "\n",
    "Stationarity is an important part of ARIMA analysis. There are several avenues we can explore in order to attempt to make a time series stationary, but one of the more common ones is differencing. `pyramid` provides a differencing method that can be used for this (`diff`), and also provides a function that can estimate the number of diffs required to reach stationarity (`ndiffs`). Before we can difference the actual timeseries, however, we need to address the pretty apparent seasons in the data. We can estimate the order of seasonal differencing using `nsdiffs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pmdarima.arima.utils import nsdiffs\n",
    "\n",
    "# Test for best nsdiffs (assuming m is 3, as previously shown)\n",
    "nsdiffs(y, m=3, max_D=5, test='ch')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`nsdiffs` shows that we do not need to perform any seasonal differencing, so we can turn our attention to differencing with `ndiffs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pmdarima.arima.utils import ndiffs, diff\n",
    "\n",
    "# Test for best ndiffs at the p < 0.05 level\n",
    "ndiffs(y, alpha=0.05, test='kpss', max_d=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ndiffs` asserts that we should only difference once. Let's look at how the time series looks after differencing once:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEICAYAAACzliQjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wcmMxLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvjNbHMQAAIABJREFUeJztnXeYFeX1x7/nbmXpdV06ShEMAoqKAZQuoImxa0xi15goSWw/jD0iwRq70diNYosdBKQoKCgsVdrC0vuCdGFhy/v7Y8qdO3fm3pl7p957Ps+zz95p75x5y5kz533P+5IQAgzDMEx2EPFbAIZhGMY7WOkzDMNkEaz0GYZhsghW+gzDMFkEK32GYZgsgpU+wzBMFsFKn/EdInqdiMbIv/sTUZnmWBciWkREB4hoFBHVIaLPiWgfEX3gkXztiUgQUa4X92MYN2GlzwQKIcQsIUQXza47AMwQQtQXQjwN4EIAxQCaCiEuspouEZUQ0WdEtFVW4O0dFdxhiGggEc2QX27r/ZaHyRxY6TNBpx2AZbrtVUKIapvp1AKYBOACpwRzmZ8BvArgdr8FYTILVvqM5xBRLyJaILts3gNQqDk2gIg2y7+nAxgI4FkiOkhE4wHcC+ASefsaq/cUQuwQQjwPYJ4D8l9FRCtk+dcS0Q2643cQ0Tb5q+Ja+cuio517CCHmCiHeArA2XXkZRgsrfcZTiCgfwCcA3gLQBMAHMLG+hRCDAMwCcJMQop4Q4jIAYwG8J2+/QkT9iGhvgr9+LjxGBYBzADQAcBWAfxHRSfLzDQdwC4AhADoCGKC9kIhGJ5LXBVkZJgbumGK8pg+APABPCmnipw+J6JZUExNCfAugkVPCWbznBM3mN0Q0BUB/AAsAXAzgNSHEMgAgovsBXK65dhyAcd5JyzCxsKXPeE1LAFtE7Ex/G/wSJhWIaAQRfU9Eu2XrfCSAZvLhlgA2aU7fFJcAw/gIK33Ga7YBaEVEpNnXNtXE5CGeBxP89U9f5Jj7FQD4H4DHABQLIRoBmAhAeZ5tAFprLmmju/7vieR1UlaGMYKVPuM1cwBUAxhFRHlEdD6AU1NNTB7iWS/B3yzlXCIqBFAgbxbI28qx+4noawu3zJfT2AmgmohGABimOf4+gKuIqCsRFQG4Ryfv2ETyauSJyPLlSZtUKPeHMExasNJnPEUIcRTA+QCuBLAbwCUAPvLo9ocBKNb0SnlboQ2A75IlIIQ4AGAUJOW+B8BvAXymOf4lgKcBzABQDuB7+dARm7KeIcs3EdKX0GEAU2ymwTBxEC+iwjAAES0CMFgI8ZPD6XYFsBRAQQqxBQzjOKz0GcZhiOg8SBZ6EYA3ANQKIX7jr1QMI2HZvUNErxJRBREt1ey7n4i2yHOjLCKike6IyTCh4gZIY/nXAKgBcKO/4jBMFMuWPhGdAckf+qYQ4hfyvvsBHBRCPOaahAzDMIxjWLb0hRAzIXW8MQzDMCHFiYjcm4joDwBKAdwqhNhjdBIRXQ/gegCoW7fuyccff7wDtw4ne34+is17D6N14zpoXMSj8BiGscb8+fN3CSGap5OGrY5ceTraLzTunWIAuwAIAA8CKBFCXJ0snd69e4vS0tJU5M0IHp9Shmeml+OWoZ0xanAnv8VhGCYkENF8IUTvdNJIa5y+PHNhjRCiFsB/kEaQDcMwDOM+aSl9IirRbJ4HaTwykwQeJcswjF9Y9unLc5kPANBMnu/8PgADiKgnJPfOekhD1RiGYZiAYlnpy3OZ63nFQVmyhpipxhiGYTyE595hGIbJIljp+wD79BmG8QtW+gzDMFkEK30fYJ8+wzB+wUrfB9i9wzCMX7DS9xE2+BmG8RpW+j7CBj/DMF7DSp9hGCaLYKXvI+zeYRjGa1jpMwzDZBGs9H2EffoMw3gNK32GYZgsgpW+j7BPn2EYr2GlzzAMk0Ww0mcYhskiWOkzDMNkEaz0GYZhsghW+j4geLAmwzA+wUqfYRgmi2Cl7wPEgzUZhvEJVvoMwzBZBCt9H2CfPsMwfsFK30d42USGYbyGlb6P8LKJDMN4DSt9hmGYLIKVvo+we4dhGK9hpc8wDJNFsNJnGIbJIljpMwzDZBGs9BmGYVJECIHDR2v8FsMWrPQZhmFS5Klpq9H13knYd6jKb1Esw0qfYRgmRT5ZuAUAsPvQUZ8lsQ4rfYZhmDQRIYq0ZKXPMAyTRbDS94EQGQUMw1iAQhRpyUqfYRgmTTLSvUNErxJRBREt1exrQkRfEdFq+X9jd8TMLEJkFDAMk2HYsfRfBzBct280gGlCiE4ApsnbTBJCZBQwDJNhWFb6QoiZAHbrdp8L4A359xsAfuOQXFlBmPyADMPEE8Y2nK5Pv1gIsU3+vR1AcZrpZRVh8gMyTDbz1NTVaD96Ao5Uhyv61gjHOnKFpMFMtRgRXU9EpURUunPnTqduG0pCaBwwTFbz2ux1AIBDR1jp7yCiEgCQ/1eYnSiEeEkI0VsI0bt58+Zp3jbcsIHPMJlFmJp0ukr/MwBXyL+vAPBpmullFWH0BzIME27sDNkcD2AOgC5EtJmIrgEwDsBQIloNYIi8zViEffoMw3hNrtUThRCXmRwa7JAsDMMwgcTMPgvjtzpH5PoIu3cYJtyE8VudlT7DMEwSMsk+Y6XvI+zTZ5hwE8Z3ASt9hmGYLIKVvo+wT59hwkEmfZSz0mcYhskiWOkzDMMkIZM+ylnpMwzDJCGZeydM7h9W+gzDMFkEK30fCJFRwDAMkrt3wuT+YaXPMAyTBHbvMGkRIqOAYZgMg5U+w/jIY5PL0HfcdL/FYLIIy7NsMs4Roi9BxmWenVHutwhMOoTws50tfYZhmCSEqaM2Gaz0GYZhkhCmjtpksNJnGIZJm/C8FVjpMwzDZBGs9BmGYSySCb59VvoMwzAW0fv2w/gOYKXPMAyTIuHx5EdhpZ8GFfsr0f3+yVixbb/h8QtfmI3/zFyLGWUVaD96Ak64dxIOVFZ5LGWsPC/PWuvb/Rngzo+W4I4PF8ftH/nULB+kYYLC0i37cOL9k7Hr4BHX78VKPw2mrazAgcpqvDF7veHx0g178NDEFRjzxXIAwM9Ha7Bsq/ELwgtKN+zBmAkrfLs/A4yfuwnvl26O27/cxHBggoXesnfKvfPizLXYX1mN78p3OZSiOaz0GYZhfMbLvgFW+g6QSYEbDMOYI1xu7F7oElb6DMMwaZKusvZyKCgrfQewW2D8ZcAwjF+w0ncAO0o8jON6GYaRcNteEx4MAmWl7wOZENXHMAxADjVm7sjNcNi9wzCZgdMdu9yRm4FoLQO2+MPHz0eqcc4zs7Bs6z5f5fh+7U+44IXZqKqp9VWObMOKUp6zxn7ZOPXFYAVW+g6QanmxxR8+5q7fjaVb9uPhSWW+ynH7h4sxf8MebN172Fc5sh0jZa2Uzba9lbbTY0s/JNgtKLbww4tSdG6P106GcnvioQGeYqWjNejtm5W+D7CFH168/AxPhKr0gyEO4xBelCcrfR/hBhs+opa+r2KocB0KBkbVIZXhl+zeCTjptregKA7GOoqS9WI8dSIU91JQvjyyBgvFnorLjYdshgSrzV57HrfRcKM0aL9f2LWqT58JKqnUES+qVa4TiRDRegAHANQAqBZC9HYi3UyHXwBMqihfGlyHgkdKZeJhOTqi9GUGCiHcnww6QHB7yz5U944LJpkQwrK7hkfv+EMmeGTZvWODl2etxY3/nY+7Pv7RlfT3Ha7CNa/Pw/wNe3Ddm6WorKpJes23q3dZkue+T5dakuGuj3/ErNU7LZ2bjagduXLzn7R0GzrdNRFvzVkfd+79ny3DjLIKfLxwM56cuspRORTlE0mi8z9ZuAVPTl2FOz/6EbPXRG2yCUu24ZFJKx2VKVt47bt1eP27dQCi9eHaN0qxRRczoX9BHDxSjWten4dt+w7jvXkb8fzX5QCAdbt+xkcLtkjXeOA3dMrSFwCmEJEA8KIQ4iX9CUR0PYDrAaBt27YO3dZbtKtOPXRed8fTf3/eJkxbWYFpKysASAp9SLfihNf87pUfLMnzxpwNlmR4+4eNePuHjVg/7mxL52cdOkv/j/9dAAC459Nl+P3p7WNOfX32eryuWVXtr0M6J0xaCOuuASGsOfX/+t4i9ff4udFy/fM7ktx3DD/e2g2zHK0yfuBzaSW8K/t2UBX7xt2H8MiklXjq0l6mRTJhyVZMW1mBpl+tUldP+9OAjrj7E3eMSDOcsvT7CSFOAjACwJ+J6Az9CUKIl4QQvYUQvZs3b+7QbTOLWr97B5mkqB25LqRtJ0127/hDoiaaaklo0/RCAzii9IUQW+T/FQA+BnCqE+lmIokqRg0r/cATlI5Tof7nOuMFZn0tiaqDVVeN180+baVPRHWJqL7yG8AwANYcyEwM+sIPioJhoqhF4lJHrhvnMs5hbRoGew03Jk0PitUJn34xgI/lB80F8I4QYpID6YYGp6yt2lpuyEFHadBuWNi23DupXMQEEq/f32krfSHEWgA9HJAl62H3TvBxc8imHZT7c43xBitfVnoLP5XgTS/gIZs+YGYl6g19du8Ej+iQTeex8yJh944/uJLtMd4dXi4x0DiulLkhhwY3lK6dBq925HKV8QQrfnr9Gfpts7LyetQeK30NL89ai+/KrQcVK2X1+eJtuOGtUrw3byNu/O98XP7y9zH++bU7f1Z/f7Z4qzrM7pFJZWg/egJueW8RnvhqlaF7RwiBhyetxPKt+9V95RUHcd+nS3HVa3NNZfvyx214f94my8+SySzdsg+PTU6+6MmEJdvwfukmrN15EN3vm4wvf9ymHhNC4NHJK7Fi+wEAwP7KarQfPSHm+klLt+PduRtN03/l23Ux2/rguzs/+jFpQN6EJdvw4fzNcZ8aL89ai4cmLMebc9YnvB4wfmEdPFKN0f9bgoNHqpNeb5UZKysMg9a84rvyXeh+/2Qs3rQXy7fux6OTV6rPvnDjHjw1dbXltJTr3i+NtqmV2/djdcXBmPOmLNuOdbuk9m5Flb/zw0aUbtijuU/8OVOX78B/v7cWZ2MFJ6dhCD1K8JXdwKTDVTWYvGwHJi/boe5bsX2/4blvztmAP555XMy+jxZK0Xg3DojdTyAcOlqDF75eg7fmbMDSB84CAFzx6tz46D9dCP+Nb0vBNxef0sbWs2Qiv3nuO1TXCvxtaGfkJAhhVQKWjmteFweOVOPGtxeodeFoTS2em7FGPbdc19gB4I//nQ8AuPRU4+DDB79YHrP9QWnsS/mjBVvQq21j/L5Pu6Qy1s3PARD9OtAGDv5BFySmZ++hqrh9r367Du/O24TiBoX429DEQWRWuer1eQAQF7TmFZe/LAUunvvcdyjMi6CyqhY3D+qEwrwcnPf8bADAX4Z0spXmk5oXxfAnZ8Udv/6t+bbS+7uFaPpr3ywFAPwuQb2wA1v6LlFdY/6eN/fpx+832ldjMMqHP/PNsdtBXu1R/hqXo8Wx3WncN1E9y9RqVCsvV5uqS9bJKaz9Dqhjpe8S1bX2F6w2GrKpTqGrqSdGL41MbaxOYlWheuVjNVx4w+Kt1dE7AZ2+N9OwVHd0utzskkR9N6GJyM1GJHeK+fGqBJa+2Zs+TudT9EUQ0dzMLZ2U6SNCrD5dCu/rlDDKbuuWfupllfCllqF1gCOXo7DSTxEhErePVNw7Rp/7SgPV+qKNLcT0K3WGtneVdJ7PlamULe4zvDadcfoGF/ntcsgErOZhovN4ucQQU5WKe8fQ1yr91/Y/GlqItu+WPUQVZOruHbsuHysvYaNzLLt3bEkTi1Hgd7ZYwt4aNsHMU1b6KSKQuFMokaVveo2uNRKiiiG2I8mesvCiczCTMFLwdvMmVeVi3dIXMf/t3YNL2i5erUXMwVkBJlljq66xb+kbpamMPEm2WEbCziHLnYOZrQzsdpJqsW3pWzjH8OXiwcyMGV7MhqT7zNamYXD2nm7hi9Lfvq8Sew8djdsvhMATU8qwbd9hg6vSZ8HGPbjytbn4ark0nv5IdQ3GTlxhGJDy9g8bcNa/Zqp+9vkbdscc37D7EMZONF95SBknbwe9T/+9eZswd51035wkHblCAPPW78aQJ75BmRxApLBw096YbSUAbfKy7ZiybDtGjV8Yt5DDrNU78emiLUllFkLgqamr8dyMcsMVt/4zc22cPG5SXnEA//5mTcJzlm3dh2enr8bYiStwtDr+5Vxx4EjcPrsNOJmSGPHULHW1JC1jJqxAn7HT0HfcdIwavxCl63cbXB0bkbvvcOy4+5dnrcWKbfsx3iBQ7L5Pl6rj1wGoQT+z1/wEAPhO/q+ntlbgkUkrsdMgb4z42cEgr2RUVkntWLmnEAJPTzMOvHpuRnlMAJ0QAle/Pg9jJ67AuC9X4rkZ5ViwcQ+OVtfi9H9OQ7+Hp2PPz/G6yogP52+O2d518CgenrQSNbUCX5dVYPRHUhubtnKH0eUAgH9OXIkj1TWuthtfgrN2HjyChyaswKMXxc7TtnTLfjw9vRzfr92N9/94uuP3PV8OyPi6bCfWjzsb787dhJdmrgURcOeIrjHn3vWxNDv0zNU7MbBLC1zwwpyY44Mf/8Zx+fQfB18u3Y4vl24HEPt5aaZOLvq3JOOvnv02Zv8FL8yO2b55/EIsuGcobtAFkvxNs7LT71+Ron3P7dkqocybdh/GvzRLAeoD2x6auAKPTinDqjEjEqbjFBe8MAf7Dlfhqr7tUZCbE3NM0cNnPx3Nn3ZNi3D5acmDXux+BSU7e8U24+A9ANi+vxIAsGXvYSzdsg/TbxuQUJ5xX8YaH2MmrACR8YtKv4La3Z8sxe/6tFONi/ma6FAt36/9Cc9/vQZl2w/glStPMZVd4dkZ5UnPcYrxczfipZlrESHC6BHHY+u+SjzxlfHylM9Mj5Vr3+EqTF9ZgenyanUKY8/rjm37pHJ44PNlKcl118c/Yu2un3Fahya48rV56v5dB81fIgePVOO9eZtcbTe+uXeqDNwfSufn0RRcI6mgWHmJ/O81KfjmUyWRYolErJ0HwNB61bLbxHJJZWbnRG4PRc5k8jjJ4aPm0xgYucCMRkwZXmvb0rd3vhlGXx1AbFkZ5a/TroUq+YZW22aVh2WutF/FpWpninKz8tfG2VSlOOW5MqWGXdegMtzbrXYTKJ++oiSS+a8du5+sBBLdzsvJkBLdKZLE0ndCzFQ6kRJd4cfyAE70bRhh36fvzcML4c29atW2aa1x+unOtjO6zYqcqaojNbDSZgpu9635pvSNesOjwxM96ilXCiXB7bysvIkUS7I8CeKIDD/X/HX61l6N3gkqXhtk6WCnLZjVUe1jpqqP1LRtXu523QmUpW8UfeomSt4mGo7l5YiWRJYxJRun74SYKaSRKH/8UPoJh66mkW5QF60XEJ5YJoq3w7Kl72N22bq3hXNTVUe1qel81w04/yx9g33K8ESvFg+x8iJ2pfKapJnQpx8zeifxOPJU88/pRw1aw0/rBe6TTz8o1KptM/imvp2sNzW0NM+ZqhFqHGOThkwOEShLX2koiaa/dfR+SK71vfRLJ7pVjoXROwqp+yCd9un76N5xeFI6u/XAW5+++xhN/BdUjCOdjXPJSh1Ntz3ZtvQzVukb5ITdzqJ0iVr6Cdw7bjQpk9slskSTZYn22lStMadfcH6u82441UBahr7NjtwMs/SVV0s4fPrWsaT0U7X0U7rKfWMpEIuoTF2+A9v2HcY9n0rjYatra/HMtNW47oxjUZiXg7nrduNAZRUGdy22la4QAn//eCnaNKmD6/ofG3Ns0+5DamDKh/M348KTW6vH7vwoGqj07txNaNmoTqqPZohZsEei6XpWVxzEY5PL0LZJkWFtGvhYNG5AOwxtsS4wS0G7KpTCJS/Oidt32weLMWpQJzStl4/XvluHGwd0RE6E8MDny1CvIDfuBX2gsgpvztmAtk2K0KgoLy696ppa/OXdRejXqRkuO7UtNu85hOkrK5Iu/GHGa9+tw9Buxfhk4RZc1beDJmhJ4IWv12D9ruiqZVOWbUeHZnVjrp+ybAdOaNkAJ7drYpj+W3PWY92uQ9i+33rA4JAnvsG1/TrYfRRDDh6pxjWvz8MvWjVEXk688tmw+xBmr7G+2psRf3g1dgW279f+hD7HNsXmPYdw5Wvz8MbVp8YMspi0dBuWbtmPMzo3xyntG+PlWevwm16tMHnZdggAO/ZV4ptV0XHvK7fvx/Kt+3H+Sa2RDpVVUtDSDWceh/zcqL36+uz16u/Hp5Rh8eZ9cdeaGSD/NAmwVAI4AeB/CzYbnpMMZbEaff4mwyywbPOeQynJoYf8CL0vKOkkbnrqQzx+sRScpV92TmHU4E64ZWhn9bjdFa1K1+/GhXLA0j/P7x6jzNs0qYNNu6MNOTdChotneMlZJxTHrL5lRr2CXEeXtUtGcYMCjPhFCV6fvR5PXdoT5/ZsZVpmvz2tLd75IT4SVCm79+dtwh3/W6LuG/jY11i362csuncoGhXl25Jry97D6Dtuurp95S/b463vN6CmVmDxfcPQ44EpltNaP+5s02fKRtaPOxvd75+MA5XVaFG/APf+qhtuemchzu5eggkag2HiqP4Y+fQsnNK+MeatNw7s0qaZDs/NKMejk8tw99ldca1sxO3YX4nTxk4DAFzbrwNe1i1JGXa09fLYZnUx4/aB84UQvdNJM1A+fT3J1gtNxhFNcMMhXdDO3p9jQ9f9VvhAcOfI2XuoSn3JHEkSMJIs/P7Q0djjynQcVoOktOjnN6qsqlHzxij4j7HHgUqprPYeqjL16R+pltqVfioIN1DqllYvBKHdesVeh/I40Eo/XbS+MTtRen5hVUS/nySdl47+GRV/aSopxqcV/c1K30FIO04/VusrZeClHUJJRrIxifExOMv9e2jrQ1DHWWvxYoZFJ3By5KNSDVJJU59fWmXg5dQPmQ5BO2Qz9piy34v2lewOdtdCzlayxtIPQ4WwKqHf0bfpNPB4RS3vT+GZ4ix9ze9kbijGOtrJ2+IsfbkQ/KqR2urEX3fWCFRwltNoK2IIdL5lZer3s6TjKYuXnUz2W0otZitCpO5hS99ZzHz6ijHlRZ1Mdo+q6hA08gCQ0Za+1qpMpaPQa/xW5lZJx9LXX6sokVTKJ5FPny195xBCG2hEcccAb92nZlOSsKVvDfbpB4iwWPpu+PRTUfp6OUizz66lH4aOfr8QMJ9wrdZLS99oamzNjY+w0reEb8FZ75duRufi+gnHZr80cy2mroiOW/9k4RYUNyjEsq37MHP1LjxzaS80LMrDmp0HsabiIIadcIx6bnnFQUzRjHlfuS12FZoDHo5zt4qyolUyvFpvQOFIda26KtAdHy7Bt6vN5fx00VbD/Q98vgwlDQtjFvwY+sQ36nzxl770PU47tglObtcYv+rREl8t24H+nZvh0JEa3PvZMpzRqRkuOKk1ZpRVYOHGvdh18AhmropdqWvL3mjchd08eiHJalvZxoafokFtR6tr8X//k2JcPtCtDvX8DCnfNu62Fzi0ZudBvF+6CYC0ZsXmPYdRkBfBHcOPR9O6+Xhm+mp0aFYPrRrVwbKt+9Ckbj5e/GYtAOCRSWU4Wl2LBRv3YosmYGnCkviAw7Bz6kNTHU/Tt+CskiueTDudYd2K8dIfehsGb3GgTXjp17EZvi3fhR6tG8ZEV/bt2BTflRsv56fnP3/ojeveLHVLxIxnSNfiGIPLCdb9c6Q6wsqsfRY3KMAlvdvg6enerbwVFprUzcfCe4dldnBWMn6yuHYlEy6U5QL167Ful5evs4ISNMSkxta9zq9TbcWDtmP/Eew55H6gVxhxykAPtdJnMptIxLjT0Ao8eic93Pj+D0O/WpBxqtsp1EqfK1Fmo59i206sBSv99HDD7cvtNT3Y0gePuMhUlMqdoxviZWeEDw/fCx6W55byfaKRYOJUroRa6Ychypaxjxr9GTGO/rQCj9MPHmzpp4dT2eeI0iei4URURkTlRDTaiTStwMZcZlJjZunbce9w5UgLN/Qz6/z0CIx7h4hyADwHYASAbgAuI6Ju6aZrBb3lx+6ezEBdQU1v6dso3iNVrPTTwQ0XC1v66eGUenMiOOtUAOVCiLUAQETvAjgXwHIH0k5I2Y4DMavMPDl1FXJzQu2xYgB1cZsV2/bH7NcP4UzEUyarDzHWWLXjoONpPjV1NRrUiV9NTc/3a3c7fu9MwKkXsRNKvxWATZrtzQBO059ERNcDuB4A8o/p6MBtJZ74apX6mwM6GCa4WF3VqrzC+RdOJhAon74VhBAvCSF6pxtNpmfN2JE4u3sJAOCJi3tgzdiRWDN2JP51SQ8nb8P4RN+OTT25z+1ndUGdvBxP7pWtXNy7tdo+v75tgOl5a8aOxLvX90mYVuvGzq5bHQacUvpOWPpbALTRbLeW93mCdix3bk5E3c5jN09G4FU5FuRGPJkEMJvJ07TPgjzzcs2JkOEi8FoKs/AF7ZR7x4kWNQ9AJyLqQET5AC4F8JkD6Vqm1mS0BxN+9It2uAUbCe6TqzHQkuU3JSn3bCyvwFj6QohqIroJwGQAOQBeFUIsS1syGyhBO9p6oJ/3mwknXpViToR4SKHL5ESiDTSZ0k5mwOXq53jOApwa/eTI1MpCiIkAJjqRViqoQ/zY0mdSRFp1i7W+m+RqXDb5yZR+EqWehTqfI3K1RC39LKwJjCNw1XGfnBj3TjKlnvh4MvdPJhK60TtuUmMSts8wVuG64z5al0wyA40tfffICKWvROJyRy6TKtqlFhl30CryZJZ6sn5aduWmTqiVfo82jQAAF/VuDQDockx99diJrRv6IhOTPhed3Bp18nLQv1MzXHBya3X/2SeWoGld8+U106Fji3rs0XeZXm0bWz43mVJnpZ86vq2RCwAv/f5kDDvhGHXptJUPDsfx90wCELv04b2fLsWbczbgrpFdcd0Zx8alc27PVji3Z6uYfW2aFGH9uLPVtG8d2hkTftyGldultXLHX9cHnYvr4eQx0hqUa8eOxLF/960v2hLnn9QKHy3wLATCkOIGBZh5x0D8ZfwiTFq2PeZYnbwcHK6qQa+2jbCm4iD2V8auQ/zj/cMQIcIJ9002TV8p90cv6hG3T6G2VqBGCOTlROKW3RvFO9BbAAAZAUlEQVQ1uBOenrYafx3SCU9ONZ6K4cpftseNA47DaWOnAQDKxgxHQW7wx31r24ebTBjVD3k5EQz710x1n7YtpYJRHuvLVYvyVdC0bj7m3zMUa3YexODHvwEgLRuYrs7/7Wlt8c4PGw2PHdOgUF29LQi0alQnZv3ndAmUpW/29k7W028pbZ0TMEKxn5hh8OkGYRhqbiRiqiDrFkg2RGFujmF+5kTIkQCoSIRMh/wJiyO5GmrmgAmDwgfgWfAYgRwfFGE3j/Xll2zbLokeL9NHcflq6etxU++STsnXCu/GgDtFEL5oFWVgNGZY0cMRk3e0F5/k0eG7ic8zfGkEvK179dIn8r9/Q//SId3vILSFsBIKS9+NtGuFCF3FCYK4SmM00glKR7pZOXoxpFaZfjZRR6EQIpTDe72tr/5qfX35aKUhSl9XJHqB+v3Cc5tAWfpuVuocii3mmloRCHeJHYLwklKVvkHLUFw6RMY5m0ME4fIzKJZ+orwKa5v2qviDYOkbGWlagtAWwkqgLH0z68yJ+ie5d6LbNSJ8/p0gvKQUa95IKURUS9/42ohDPv2EKDEbCS19s0uD/TrwMiDJ75zQP2qskUFZNXrH6UXqA6X03UTfiVhbG0L3TgDkjSTw6SvKPlGDdPsRVEs/hWv9tm6T4ZmlHwDjQl8W+lWjQuidCwyhUPpOlK/eRyi5d8JFEJS+0v9ppB+VF4J+ZFTMOS4/hLBi6ftuxwafoLwAlWLUGhmO+PQz0P1nFV+Vfufi+slPAnBG5+YAgN7trQd36DmpbWOc2yM6lr9zcf0Qzt/hn7xDuhYDAM7rJQVLjZQXrtFyTb8OpscUlCz/w+nt8Ps+7RyV8bJT26K/XFdO6dDE9LwBnVuov08/1ptFWpzAbnW9uHfr5CcZ0KJ+AVrULzA8dtYJxejdLr4d9uvYzDS9y09ra1uG+oVSd+Mlp7SRZSpUj116SpuEdcwKJQ3DswiL0y8hXzty2zera+m8Mzo3TzmAZvVDI1BTK1CYl4OOLerhoYkr1Hv/fKTa8JrHL+qBWz9YrG6vGjMCne/+0va9vaJF/QJUGKwf+8TFPfDo5DJs2xcfaPLMZb1w8/iF6N+pGW4b1gXnPvdd3DkrHxyOCElusRwiVFbXqKtLXdy7Dc7t2RJd7o4GC11+WjtceHJrFOTm4KEJK2LSWv3QCADSF8CqMSPUeVjuOacbAOlLLJWpY8sfGgEBoLpGoCA3gkhESj8/N2rPXPnL9nh99np1e0i3YlUm7dQdVu5emBdBZZJF16fdeibaNC5KWmeW3D8MJ94/JeE52rqnNVKUIDQA+PCPp+PCf8+Ju/bhC07EMQ0KEy4jevtZXfDo5LKYfY3lqGd9vS9/aIQ8G2l0kkNFpBwiNbhx1ZgRamdwVU0tivLtt9vCvBysGjNCnZitSd18rBoj1aG8HAIRoVNxPfz62Wi9LRszHLmRCI7TBVku/8dZeG5GOZ6bsQY3nHksbh3aBZ8v3goAOOfEEnyxZFtCWcrGDFfree92jVG6YY/heaV3D8F78zap+bn8H2fh9g+XYEKS9L0mUKN3EpFqAE1eTgRmi+yYWU76VXm0CsRfjNWSWaNKJLfyjAW5ETWoyuyc6H1izzMqE7Ny0o6L18qVr52EK4UvmVw5Xa2o+uc2m9ExlYU46hfmobIq8QLtRfk5lupMoYU6bZZOvuaZck2eg4hQZFK2Colmu9TfW3ufRENeY8o3jbajvza+XGO3zepeUX6uem5BTgT5uRF1mmcrX/u5msCTJgmmAWlWL/brqCg/FwWB0R1RgieRi+jL16zDyqnFCrwiUcU1OxJ9Rg9G1PiM1eK0MkrCSlZZ7QhNJ9+tuiZDVpXdR843RZHX1MZ/tenzTJvTfowacroMs0vp6xqjUn5xw8M8kscuqRS+6TBYtcPTnYocpBeJk+Vp5bmsjiwJUBaFEjvtQTlXyXPlS6WqxsKLXlNQZtHmYSIDHsE6Zg1Wv9vpcbFuY6Y8Ej2Gdo6aTFc+li19C+dYyi2rSj+NN2OQXqp+YWcUllLflXxT+pRq9GNBDYiZoytJxodBd2SX0tdvq5a+LuQ7oOVmJleiemh2LDpdQeYrECeHaFqz9C26d9KRI+Nf1e6g5Jvi06+2oPRjrk+q9FOTy0uyS+nrCkypAHGWfmAdPMYk9OmbHFKe0T0fZXCUkpMN0ZpP3320xZatLw9b7h3dtuLTr65JPBJLTyYEhWWX0tdtK4pPr/cM+nYCgZMvI9XAoXA3fCs4+cltxSVj2dIPQLaHzcBJFW0fFpDM0jfPk2Sr87mRm06XkW9Kf/gJx6i/TzYI9nADfXkpb/vzeklBWyUNpQAQ7QpcfmLVqjhTDkgywkyhd2pRT7q2U3NHlI+SnsKve7QEAPTweQWz+gW56GMxAMvo3VAvyZBHI5T8HNn9mCTnpeHT1/w+pmGh6XnJaN24KGb7V3K5pUKvto1SvjYVWjQwDiArNtivN/DaNJGee9DxLWLOG3x8C5xzonkeeB3QOVgnnxP4Mk6/W8sGeOa3vdTt8df1wVGbn1mpoC+wnAhh8X3DUFce5z7jtgGoFQJF+blYcM9QnPTgV3FpnNG5OWau2omx53XHr3u2xCUvzsGyrfvxyIUn4o4PlySVoUfrhli8eZ8leesV5GLWHYPQ4x9SAI9eKV14cmv83/DjcaCyCq98uy7ueoFoJf/4T7/Eec/PVo91LWmAhfcMReO6+di0+1DMdYOOb4FnNeVjhS9G9YvZvuvsrvjL4E6ok5+DI9U1ttJyihX/GA4iKd5g4T1DUVSQY2m0hpbSu4egulZgwKMzsOvgUTU/n7q0J/7y7iLDa5R69tSlvTDughqc+cgM7DlUZXhur7aNsHDjXlsySfeQ/l92ahsUNyjE/LuH4HBVDS558Xts2XsY71x7Wsz555/UCvee001tZ6c+JK0a1rpxHSy6dyh6/kOq6/fKwXKp8N71p3vSjhVa1C/EonuHojAvJ6ZD9pvbB6KyqkZ9JiC+7bRqJD13wzp5GPflSgDSym6FeTmIEMUE82kxM8QK8yIx97m6b4eEslsJzFPk6ffw9ITn2cUXpZ9DFBes41cAlHYFJW0wUuOiPKPT0UwOzsjPjaBeQS4ayee1tBjW3djmGq8NTeQApFD15vULcKDSWKEAUYuwUVH8fc1kaVgnLy4QKxn6wJicCKmy+1W2dTRBa8qz2jXclTpRvzAvRuk3qGNeLso5UmBgBEX5uaZKv36heTqJUL7glHxvKgcGNatfgC17D6NQF7DXvH6BYR0ApLrRuCgPew5VpbXOgB/t2OiZCvNy4gILVW+mxvDTX2ulLMxcd0p7Ub4o6hUkDrxrYOFeijw8Tt8jkn3GedNZl0yG5FGFyrFEfu0g+JbDguouS9AQvRgCnKzMrBap3v0XhiGHqeDUY5mN01eCHaMxj8FtVKz0PSadqpBKvVXul5lN2Tv047wTETcE2A2BzDBfLCAh4Zt80B5mgzbskizYUT0vvdu4Cit9m6TbgJ1sXHaSSmTpxMcp8CvCjOhLNMGXk27bSpCcUxi5MRioGZPuSDUz75dq6aeVujFOp8lK3yGsDquyU+XiVw+yl5aIWR3MupJi4lFyL6K6y+xc64Z7x6b70SwaXfFWZcmLPt13oZlPv1Y39DPI71xW+jaJi+q1qTK9rgyqZZodbdp9LLl33BfD7F5czsZowlLSwlTpKzcIQQGw0rdJ+kWaerUzsxgTd9Imv19YJpzzE/2EXQndZXGdo27IY1IXUvRdZ7o7qLY2tXzRY670dZa+g9/PWTF6p6Mu0MdpEq3ykwxl9a5jm0sLwJx+nBT407KR8yvxDOoiBWaUNCxEfk4EvdvFrgaVrAJ3aFYXA+TArURDP5UApBNaNgAAw5WRzAjifOGpMlReXMWIgV2kfBwmBxW2aVJkGAQERKM9FQZ3TS3Apm/HaFBZU4tDfft3kuRsbrLylRlKEJB+qCNgPnw5DCiL/vRoIwWOdS1pYHhefYtjec18+krd0Q/e0bdZK7RvWpT8pDQI3CIqpXcPSWmlHavMvWuwpTGyZvz21LY4o1NzNaLvxjOPw7k9W8ZFNgLAt/83EEeqa/G/+Zvx/NdrAJgr6ll3DET/R2YAAL6+bQCO1tSijZzmpzf1RXWNQEnDQvTv1Az/mbUWb87ZoF6rWGltmtTBe9efjtwIYX9lNTq2qIduJQ1wdb8OMcvNLbhnaMy96xfmYc6dg9CsXgG276tE68bWX2Cldw9Btc2Ap6Dy7G97xawENuIX0Yjae87phhvOPA4lDQtx+Wlt0bpxEabfOgCVVTU4ecxU9bzbz+oSt7jHP879BUYN7oTT/xkNsvnm9gFJ5XnlilMASOWljH//bvQg1NQITFm+3fCa24Z1we/6tLO8HKBSH8ee3x23DOsct6CO9t5hY+E9Q9UX8K96tETPNo3Udmt2nnZfL4PgzIiJ1n/84h4A4r8ILzu1Dfp3aoa6BbnIzaGYgKz5dw+JqTsKE0b1R2VVfEDjp3/ui54PG97eFoFT+vrVZ5xGq/xSgYhiKk4kQoYKH4iGuGu/AsyMc22a+mUktTK3aVKEVnJ6cRHGROq9WsgGTW5OJO4rxGj1H0VJGDWKRKQaXBRE9AFm2rqozUelXOsW5MYpyeIG8fUrLycSp4TNAqUUIhS1urXl1UpXlno3Qk6E4s6xgpGM+nuHDX3woVndNgpSNAtcNDPalLqjd8Hq9YWWpia6Lr5eSWmWNEpPdymE8xUeYpxwnTrVKcUwjD38nEzPqX4CVvoeoH33Z/qMlow1vOo3TTZclOujPZIPk/ZEjLRgpe8xjlj6Qr8dgprGxOCUqs3wQTehxcnRUIEavUNE9xPRFiJaJP+NdEqwjMItpcwN3lWyQaFmwzM6SbL8ctP8cqqsnOjI/ZcQ4jEH0skKnPHpy2OOLUy4xgQTLrNwkswdpo7ecbB4eRqGEOK0T9+NisWEC6sfj+yz9xY3VyJzqiSdUPo3EdESInqViLxZAstDjndgFa32TaNDMLu3bohTO5gHbFgJdlKC1zoXS//rF0ofbInSZfxHCegD0m/ASlrJVnlLpoQahTjwym2MgqSSGVqKvujY3DzANFE7NYqRUQJACwwC51IhqXuHiKYCMFr37S4ALwB4EJIx+yCAxwFcbZLO9QCuB4C2bdumKK63zLx9IBrXlRrF93cOTjlI5YzOzTFxVH8AUqX4dY+W+OW4+NVwZo8epEYQJuKsE47BhFH90E2OLmxWrwBTbzkj6Rj70ruHxKwwxMSz4J6heHHmGrz4zVrHbeRP/txXDc6JTnQm/X/8oh649YPFltMa3LUYE0f1R9eS9IwSq0Fc2chnN/fDPnnxmyt/2R6vz14PghQ8+fhXq/D54q1x1/ymZyt0Lq6PE1oaLxM647YBhpHcpXcPweGjNYaR849f1AN/G9I5paU7jUiaihBiiJWEiOg/AL5IkM5LAF4CgN69e4dC87TVvOnTWYcUkJaIVDCz5u1M5aCvVB1bJG/8bge+ZQJN6uajxCDAygkSRYLbnTYBiK1TjPM0KMxTy0z7RdS+WV11iVU9RGSq8AFpahQjErXNwrwcR6emSXf0Tolm8zwAS9MTh2H8xwuLhH3t4SIMK2JZJd3vhUeIqCekdrIewA1pS5QF8MiNcOBlOXGVCAeZUExpKX0hxO+dEoRhsgmeBz9c6IsnzC9pHrLJMAxjkTArewVW+gwTANxWJtyHkCYZ9CnGSp9hdHjRvjPBYswmojPb2l8jOWiw0neJk9o2Mj3G7T3YKAEy2mAqK+cD0ipnTnCyjdXLGPdpJwdYtnN5VSsvCNwiKpnAFzf3ixnjz4SLYSccg/eu72M5wvnzm/ph+/5KHKistnyNajEaDBD94uZ+jisX/X1m3DbAUiAgI3HBSa3QpnEdtXyVL7W+HZvigV//wkfJ7MNK3wV+0co8OENL/cJcHKisdlkaJhVOO7Zp8pNkGtfNN11pKRWs1p90MAsSYowhIsM6cXb3lq6v6e007N5hGB/w2qfPHbmMAit9H+FmmL1w2TN+wUrfB6IR3dz0GYbxFlb6PsI6P3vRv/DZ/RJO3Jw/3y1Y6TMMw2QRrPQZxge8suvDHETEuAMrfR8QQlnjlslW9IuoMIxX8Dh9B5l1x0Bb4+65I5fR8vlN/dC0nnPj/ZngM3v0IPx08Kin92Sl7yDJlitkGAWjF3731u4HZTHBomWjOrZWzHMCdu8wTADgjz7GK1jp+wC7cRkm7IT3Lc1K30fCW22YsMBfEG4RXtONlT7D+Ijbo3d4dBCjh5U+wzBMFsFK3wfY+mIYxi9Y6fsI+1sZBa4K4aJFfWmFtEZ1whdXweP0GYZhbHLToI5o36wII7sf47cotmGlzzBZAH9VOkteTgTn9Wrttxgpwe4dHwjjdKyMO3BdYLyGlb6vsPnFMIy3sNJnGIbJIljpMwzDZBGs9P2A3biMx3BsCKPASt9HeEQFo8J1gfEIVvo+kBORWniL+gUAgLr5OX6Kw3hIka6sm9aT6kBhnjt1oEGhNCq7YZ08V9JnwgcJH777iOgAgDLPb2yfZgB2+S1EEsIgIxAOOcMgI8ByOkkYZASicrYTQjRPJyG/grPKhBC9fbq3ZYioNOhyhkFGIBxyhkFGgOV0kjDICDgrJ7t3GIZhsghW+gzDMFmEX0r/JZ/ua5cwyBkGGYFwyBkGGQGW00nCICPgoJy+dOQyDMMw/sDuHYZhmCyClT7DMEwW4anSJ6LhRFRGROVENNrLexvI0oaIZhDRciJaRkR/kfc3IaKviGi1/L+xvJ+I6GlZ9iVEdJKHsuYQ0UIi+kLe7kBEP8iyvEdE+fL+Anm7XD7e3kMZGxHRh0S0kohWENHpAc3Lv8nlvZSIxhNRYRDyk4heJaIKIlqq2Wc7/4joCvn81UR0hQcyPiqX+RIi+piIGmmO3SnLWEZEZ2n2u6oHjOTUHLuViAQRNZO3A5OX8v6b5fxcRkSPaPY7l5dCCE/+AOQAWAPgWAD5ABYD6ObV/Q3kKQFwkvy7PoBVALoBeATAaHn/aAAPy79HAvgSUsB8HwA/eCjrLQDeAfCFvP0+gEvl3/8GcKP8+08A/i3/vhTAex7K+AaAa+Xf+QAaBS0vAbQCsA5AHU0+XhmE/ARwBoCTACzV7LOVfwCaAFgr/28s/27ssozDAOTKvx/WyNhNbuMFADrIbT/HCz1gJKe8vw2AyQA2AGgWwLwcCGAqgAJ5u4Ubeel6Q9M80OkAJmu27wRwp1f3tyDfpwCGQooULpH3lUAKJAOAFwFcpjlfPc9luVoDmAZgEIAv5Mq5S9PQ1HyVK/Tp8u9c+TzyQMaGkJQp6fYHLS9bAdgkN+RcOT/PCkp+AmivUwK28g/AZQBe1OyPOc8NGXXHzgPwtvw7pn0reemVHjCSE8CHAHoAWI+o0g9MXkIyPoYYnOdoXnrp3lEanMJmeZ/vyJ/tvQD8AKBYCLFNPrQdQLH82y/5nwRwB4BaebspgL1CiGoDOVQZ5eP75PPdpgOAnQBek91QLxNRXQQsL4UQWwA8BmAjgG2Q8mc+gpefCnbzz+82djUkqxkJZPFFRiI6F8AWIcRi3aEgydkZQH/ZlfgNEZ3ihoxZ35FLRPUA/A/AX4UQ+7XHhPT69G1MKxGdA6BCCDHfLxkskgvpU/UFIUQvAD9Dckeo+J2XACD7xM+F9JJqCaAugOF+ymSVIORfIojoLgDVAN72WxY9RFQE4O8A7vVbliTkQvoK7QPgdgDvEzk/F6+XSn8LJJ+aQmt5n28QUR4khf+2EOIjefcOIiqRj5cAqJD3+yF/XwC/JqL1AN6F5OJ5CkAjIlLmTdLKocooH28I4CeXZQQkC2OzEOIHeftDSC+BIOUlAAwBsE4IsVMIUQXgI0h5HLT8VLCbf77kKxFdCeAcAJfLL6egyXgcpBf9YrkttQawgIiOCZicmwF8JCTmQvq6b+a0jF4q/XkAOskjJfIhdYx95uH9Y5DfoK8AWCGEeEJz6DMASk/9FZB8/cr+P8i9/X0A7NN8eruCEOJOIURrIUR7SPk1XQhxOYAZAC40kVGR/UL5fNetQyHEdgCbiKiLvGswgOUIUF7KbATQh4iK5PJX5AxUfmqwm3+TAQwjosbyV80weZ9rENFwSO7HXwshDulkv5SkEVAdAHQCMBc+6AEhxI9CiBZCiPZyW9oMaRDHdgQoLwF8AqkzF0TUGVLn7C44nZdOd6Ak6bgYCWmUzBoAd3l5bwNZ+kH6XF4CYJH8NxKSz3YagNWQetKbyOcTgOdk2X8E0NtjeQcgOnrnWLnQywF8gGhvf6G8XS4fP9ZD+XoCKJXz8xNIIx4Cl5cAHgCwEsBSAG9BGhHhe34CGA+pn6EKklK6JpX8g+RXL5f/rvJAxnJIfmWlDf1bc/5dsoxlAEZo9ruqB4zk1B1fj2hHbpDyMh/Af+W6uQDAIDfykqdhYBiGySKyviOXYRgmm2ClzzAMk0Ww0mcYhskiWOkzDMNkEaz0GYZhsghW+gzDMFkEK32GYZgs4v8Bp3vNzu5yK08AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_ts(x, title=\"Time series\"):\n",
    "    n_samples = x.shape[0]\n",
    "\n",
    "    plt.plot(np.arange(n_samples), x)\n",
    "    plt.axis([0, n_samples, x.min(), x.max()])\n",
    "    plt.title(title)\n",
    "    plt.show()\n",
    "    \n",
    "plot_ts(diff(y, differences=1, lag=1), title=\"diff=1, lag=1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `auto_arima` will only ever use `lag=1` while differencing, and that these estimations all happen under the covers in `auto_arima`. There is not a need to difference prior to using `auto_arima`, as it will find the optimal differencing orders for you.\n",
    "\n",
    "Let's try fitting an ARIMA using `auto_arima`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit ARIMA: order=(1, 1, 1) seasonal_order=(1, 0, 1, 3); AIC=6452.402, BIC=6484.702, Fit time=0.880 seconds\n",
      "Fit ARIMA: order=(0, 1, 0) seasonal_order=(0, 0, 0, 3); AIC=7083.345, BIC=7094.112, Fit time=0.025 seconds\n",
      "Fit ARIMA: order=(1, 1, 0) seasonal_order=(1, 0, 0, 3); AIC=6761.848, BIC=6783.381, Fit time=0.119 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/pmdenv/lib/python3.6/site-packages/statsmodels/tsa/statespace/representation.py:375: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return matrix[[slice(None)]*(matrix.ndim-1) + [0]]\n",
      "/anaconda3/envs/pmdenv/lib/python3.6/site-packages/statsmodels/tsa/statespace/representation.py:375: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return matrix[[slice(None)]*(matrix.ndim-1) + [0]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit ARIMA: order=(0, 1, 1) seasonal_order=(0, 0, 1, 3); AIC=6468.459, BIC=6489.993, Fit time=0.228 seconds\n",
      "Fit ARIMA: order=(1, 1, 1) seasonal_order=(0, 0, 1, 3); AIC=6450.426, BIC=6477.343, Fit time=0.292 seconds\n",
      "Fit ARIMA: order=(1, 1, 1) seasonal_order=(0, 0, 0, 3); AIC=6452.211, BIC=6473.744, Fit time=0.166 seconds\n",
      "Fit ARIMA: order=(1, 1, 1) seasonal_order=(0, 0, 2, 3); AIC=6452.425, BIC=6484.725, Fit time=0.507 seconds\n",
      "Fit ARIMA: order=(1, 1, 1) seasonal_order=(1, 0, 2, 3); AIC=6454.425, BIC=6492.109, Fit time=0.592 seconds\n",
      "Fit ARIMA: order=(2, 1, 1) seasonal_order=(0, 0, 1, 3); AIC=6451.633, BIC=6483.933, Fit time=0.386 seconds\n",
      "Fit ARIMA: order=(1, 1, 0) seasonal_order=(0, 0, 1, 3); AIC=6761.850, BIC=6783.384, Fit time=0.110 seconds\n",
      "Fit ARIMA: order=(1, 1, 2) seasonal_order=(0, 0, 1, 3); AIC=6456.689, BIC=6488.989, Fit time=0.552 seconds\n",
      "Fit ARIMA: order=(0, 1, 0) seasonal_order=(0, 0, 1, 3); AIC=7081.687, BIC=7097.837, Fit time=0.091 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/pmdenv/lib/python3.6/site-packages/statsmodels/tsa/statespace/representation.py:375: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return matrix[[slice(None)]*(matrix.ndim-1) + [0]]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit ARIMA: order=(2, 1, 2) seasonal_order=(0, 0, 1, 3); AIC=6454.177, BIC=6491.861, Fit time=0.912 seconds\n",
      "Total fit time: 4.864 seconds\n"
     ]
    }
   ],
   "source": [
    "from pmdarima.arima import auto_arima\n",
    "\n",
    "# We'll give it a broad range of parameters to search, though it might take a bit.\n",
    "arima = auto_arima(y,  # this is our unlagged data\n",
    "                   exogenous=None,  # if you have covariates, you can regress against them as well (optionally)\n",
    "                   start_p=1, max_p=5,\n",
    "                   start_q=1, max_q=5,\n",
    "                   start_P=1, max_P=3,\n",
    "                   start_Q=1, max_Q=3,\n",
    "                   d=1, D=0,  # we have already estimated our d and D, so no need to compute it again\n",
    "                   max_order=None,  # do not limit the order of parameters for this search\n",
    "                   stepwise=True,  # faster\n",
    "                   error_action='ignore',  # do not care if we find parameters that fail; skip them\n",
    "                   trace=True, seasonal=True, m=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Statespace Model Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>                 <td>y</td>               <th>  No. Observations:  </th>   <td>1610</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>           <td>SARIMAX(1, 1, 1)x(0, 0, 1, 3)</td> <th>  Log Likelihood     </th> <td>-3220.213</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                  <td>Fri, 02 Nov 2018</td>        <th>  AIC                </th> <td>6450.426</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                      <td>15:37:48</td>            <th>  BIC                </th> <td>6477.343</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Sample:</th>                        <td>0</td>               <th>  HQIC               </th> <td>6460.419</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                            <td> - 1610</td>            <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>              <td>opg</td>              <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>intercept</th> <td>    0.0016</td> <td>    0.007</td> <td>    0.212</td> <td> 0.832</td> <td>   -0.013</td> <td>    0.016</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L1</th>     <td>    0.1599</td> <td>    0.025</td> <td>    6.348</td> <td> 0.000</td> <td>    0.111</td> <td>    0.209</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L1</th>     <td>   -0.8480</td> <td>    0.016</td> <td>  -51.924</td> <td> 0.000</td> <td>   -0.880</td> <td>   -0.816</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.S.L3</th>   <td>    0.0557</td> <td>    0.023</td> <td>    2.400</td> <td> 0.016</td> <td>    0.010</td> <td>    0.101</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sigma2</th>    <td>    3.2036</td> <td>    0.085</td> <td>   37.733</td> <td> 0.000</td> <td>    3.037</td> <td>    3.370</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Ljung-Box (Q):</th>          <td>58.59</td> <th>  Jarque-Bera (JB):  </th> <td>276.16</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Q):</th>                <td>0.03</td>  <th>  Prob(JB):          </th>  <td>0.00</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Heteroskedasticity (H):</th> <td>0.88</td>  <th>  Skew:              </th>  <td>0.38</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(H) (two-sided):</th>    <td>0.13</td>  <th>  Kurtosis:          </th>  <td>4.88</td> \n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using the outer product of gradients (complex-step)."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                                 Statespace Model Results                                \n",
       "=========================================================================================\n",
       "Dep. Variable:                                 y   No. Observations:                 1610\n",
       "Model:             SARIMAX(1, 1, 1)x(0, 0, 1, 3)   Log Likelihood               -3220.213\n",
       "Date:                           Fri, 02 Nov 2018   AIC                           6450.426\n",
       "Time:                                   15:37:48   BIC                           6477.343\n",
       "Sample:                                        0   HQIC                          6460.419\n",
       "                                          - 1610                                         \n",
       "Covariance Type:                             opg                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "intercept      0.0016      0.007      0.212      0.832      -0.013       0.016\n",
       "ar.L1          0.1599      0.025      6.348      0.000       0.111       0.209\n",
       "ma.L1         -0.8480      0.016    -51.924      0.000      -0.880      -0.816\n",
       "ma.S.L3        0.0557      0.023      2.400      0.016       0.010       0.101\n",
       "sigma2         3.2036      0.085     37.733      0.000       3.037       3.370\n",
       "===================================================================================\n",
       "Ljung-Box (Q):                       58.59   Jarque-Bera (JB):               276.16\n",
       "Prob(Q):                              0.03   Prob(JB):                         0.00\n",
       "Heteroskedasticity (H):               0.88   Skew:                             0.38\n",
       "Prob(H) (two-sided):                  0.13   Kurtosis:                         4.88\n",
       "===================================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n",
       "\"\"\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arima.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running R's `auto.arima` found the exact same parameters with almost the same exact AIC, BIC, etc.\n",
    "\n",
    "```R\n",
    "> library(forecast)\n",
    "> df = read.csv('dummy_data.csv')\n",
    "> head(df)\n",
    "                 time occupancy\n",
    "1 2017-03-01 00:02:01         2\n",
    "2 2017-03-01 00:04:01         3\n",
    "3 2017-03-01 00:06:01         2\n",
    "4 2017-03-01 00:08:01         1\n",
    "5 2017-03-01 00:10:01         4\n",
    "6 2017-03-01 00:12:01         1\n",
    "> y = ts(df$occupancy, frequency=3)\n",
    ">\n",
    "> auto.arima(y)\n",
    "Series: y\n",
    "ARIMA(1,1,1)(0,0,1)[3]\n",
    "\n",
    "Coefficients:\n",
    "         ar1      ma1    sma1\n",
    "      0.1598  -0.8479  0.0557\n",
    "s.e.  0.0330   0.0198  0.0283\n",
    "\n",
    "sigma^2 estimated as 3.21:  log likelihood=-3220.24\n",
    "AIC=6448.47   AICc=6448.5   BIC=6470.01\n",
    "```\n",
    "\n",
    "To predict your next value, you simply call `predict`. Note that once you've discovered your parameters, you'll periodically have to refresh your model with new data as it comes in. Since `ARIMA`s use the most recent values to project into the future, you will not want to predict many values into the future; only several at a time. Here's how we can forecast the next three values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last few values: [7, 8, 7, 3, 2]\n",
      "Nex three predicted values: array([4.72557485, 4.94374601, 4.95691855])\n"
     ]
    }
   ],
   "source": [
    "print(\"Last few values: %r\" % y[-5:].tolist())\n",
    "print(\"Nex three predicted values: %r\" % arima.predict(n_periods=3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:pmdenv]",
   "language": "python",
   "name": "conda-env-pmdenv-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
